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Chapter  1 
Introduction 


The  development  of  the  Weather  Research  and  Forecasting  (WRF)  modeling  system  is  a  multi¬ 
agency  effort  intended  to  provide  a  next-generation  mesoscale  forecast  model  and  data  assim¬ 
ilation  system  that  will  advance  both  the  understanding  and  prediction  of  mesoscale  weather 
and  accelerate  the  transfer  of  research  advances  into  operations.  The  model  is  being  devel¬ 
oped  as  a  collaborative  effort  among  the  NCAR  Mesoscale  and  Microscale  Meteorology  (MMM) 
Division,  the  National  Oceanic  and  Atmospheric  Administration’s  (NOAA)  National  Centers 
for  Environmental  Prediction  (NCEP)  and  Forecast  System  Laboratory  (FSL),  the  Department 
of  Defense’s  Air  Force  Weather  Agency  (AFWA)  and  Naval  Research  Laboratory  (NRL),  the 
Center  for  Analysis  and  Prediction  of  Storms  (CAPS)  at  the  University  of  Oklahoma,  and  the 
Federal  Aviation  Administration  (FAA),  along  with  the  participation  of  a  number  of  university 
scientists. 

The  WRF  model  is  designed  to  be  a  flexible,  state-of-the-art,  portable  code  that  is  efficient 
in  a  massively  parallel  computing  environment.  A  modular  single-source  code  is  maintained 
that  can  be  conhgured  for  both  research  and  operations.  It  offers  numerous  physics  options, 
thus  tapping  into  the  experience  of  the  broad  modeling  community.  Advanced  data  assimilation 
systems  are  being  developed  and  tested  in  tandem  with  the  model.  WRF  is  maintained  and 
supported  as  a  community  model  to  facilitate  wide  use,  particularly  for  research  and  teaching, 
in  the  university  community.  It  is  suitable  for  use  in  a  broad  spectrum  of  applications  across 
scales  ranging  from  meters  to  thousands  of  kilometers.  Such  applications  include  research  and 
operational  numerical  weather  prediction  (NWP),  data  assimilation  and  parameterized-physics 
research,  downscaling  climate  simulations,  driving  air  quality  models,  atmosphere-ocean  cou¬ 
pling,  and  idealized  simulations  (e.g  boundary-layer  eddies,  convection,  baroclinic  waves).  With 
WRF  as  a  common  tool  in  the  university  and  operational  centers,  closer  ties  will  be  promoted 
between  these  communities,  and  research  advances  will  have  a  direct  path  to  operations.  These 
hallmarks  make  the  WRF  modeling  system  unique  in  the  history  of  NWP  in  the  United  States. 

The  principal  components  of  the  WRF  system  are  depicted  in  Figure  1.1.  The  WRF  Software 
Framework  (WSF)  provides  the  infrastructure  that  accommodates  multiple  dynamics  solvers, 
physics  packages  that  plug  into  the  solvers  through  a  standard  physics  interface,  programs 
for  initialization,  and  the  WRF  variational  data  assimilation  (WRF-Var)  system.  As  of  this 
writing  there  are  two  dynamics  solvers  in  the  WSF:  the  Advanced  Research  WRF  (ARW)  solver 
(originally  referred  to  as  the  Eulerian  mass  or  “em”  solver)  developed  primarily  at  NCAR, 
and  the  NMM  (Nonhydrostatic  Mesoscale  Model)  solver  developed  at  NCEP,  which  will  be 
documented  and  supported  to  the  community  by  the  Developmental  Testbed  Center  (DTC). 
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Figure  1.1:  WRF  system  components. 


While  there  are  multiple  solvers,  and  while  not  all  physics  are  available  to  both  solvers,  the 
WSF  is  common  to  all  components. 


1.1  Advanced  Research  WRF 

The  ARW  system  consists  of  the  ARW  dynamics  solver  together  with  other  components  of  the 
WRF  system  needed  to  produce  a  simulation.  Thus,  it  also  encompasses  physics  schemes,  ini¬ 
tialization  routines,  and  a  data  assimilation  package.  The  ARW  shares  the  WSF,  the  framework 
common  to  all  WRF  modeling  system  components,  including  the  NMM  solver.  Similarly,  the 
physics  packages  are  available  to  both  the  ARW  and  NMM  solvers.  In  this  vein,  it  should  be 
understood  that  the  association  of  a  component  of  the  WRF  system  with  the  ARW  does  not 
preclude  it  from  being  a  component  of  any  other  WRF  configuration.  The  following  section 
highlights  the  major  features  of  the  ARW  system.  Version  2,  and  reflects  elements  of  WRF 
Version  2,  which  was  hrst  released  in  May  2004. 

This  technical  note  will  focus  on  the  scientihc  and  algorithmic  approaches  in  the  ARW. 
Discussed  are  the  ARW  solver,  available  physics  options,  initialization  capabilities,  boundary 
conditions,  and  grid-nesting  techniques.  The  WSF  provides  the  software  infrastructure  for  all 
WRF  conhgurations  and  is  documented  separately  (Michalakes  et  ah,  1999,  2004).  The  WRF- 
Var  program,  a  component  of  the  broader  WRF  system,  has  been  adapted  from  MM5  3DVAR 
(Barker  et  ah,  2004)  and  is  encompassed  within  the  ARW.  As  a  separate  document  detailing  the 
broader  WRF-Var  system  will  be  forthcoming,  this  technical  note  will  focus  on  a  summary  of 
the  changes  and  updates  implemented  to  adapt  this  data  assimilation  capability  from  the  MM5 
to  WRF.  For  those  seeking  information  on  running  the  ARW  modeling  system,  details  on  its 
use  are  contained  in  the  ARW  User’s  Guide  (Wang  et  ah,  2004). 
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1.2  Major  Features  of  the  ARW  System 


ARW  Solver 

•  Equations:  Fully  compressible,  Euler  nonhydrostatic  with  a  run-time  hydrostatic  option 

available.  Conservative  for  scalar  variables. 

•  Prognostic  Variables:  Velocity  components  u  and  v  in  Cartesian  coordinate,  vertical  velocity 

w,  pertnrbation  potential  temperature,  perturbation  geopotential,  and  perturbation  sur¬ 
face  pressure  of  dry  air.  Optionally,  tnrbulent  kinetic  energy  and  any  nnmber  of  scalars 
such  as  water  vapor  mixing  ratio,  rain/snow  mixing  ratio,  and  clond  water/ice  mixing 
ratio. 

•  Vertical  Coordinate:  Terrain-following  hydrostatic-pressnre,  with  vertical  grid  stretching  per¬ 

mitted.  Top  of  the  model  is  a  constant  pressure  surface. 

•  Horizontal  Grid:  Arakawa  C-grid  staggering. 

•  Time  Integration:  Time-split  integration  using  a  3rd  order  Rnnge-Kntta  scheme  with  smaller 

time  step  for  acoustic  and  gravity-wave  modes. 

•  Spatial  Discretization:  2nd  to  6th  order  advection  options  in  horizontal  and  vertical. 

•  Turbulent  Mixing  and  Model  Filters:  Snb-grid  scale  tnrbulence  formulation  in  both  coordi¬ 

nate  and  physical  space.  Divergence  damping,  external-mode  hltering,  vertically  implicit 
acoustic  step  off-centering.  Explicit  hlter  option  also  available. 

•  Initial  Conditions:  Three  dimensional  for  real-data,  and  one-,  two-  and  three-dimensional 

nsing  idealized  data.  A  nnmber  of  test  cases  are  provided. 

•  Lateral  Boundary  Conditions:  Periodic,  open,  symmetric,  and  specihed  options  available. 

•  Top  Boundary  Conditions:  Gravity  wave  absorbing  (diffusion  or  Rayleigh  damping),  w  =  0 

top  bonndary  condition  at  constant  pressnre  level. 

•  Bottom  Boundary  Conditions:  Physical  or  free-slip. 

•  Earth’s  Rotation:  Fnll  Coriolis  terms  inclnded. 

•  Mapping  to  Sphere:  Three  map  projections  are  supported  for  real-data  simulation:  polar 

stereographic,  Lambert-conformal,  and  Mercator.  Curvature  terms  included. 

•  Nesting:  One-way,  two-way,  and  moving  nests. 


Model  Physics 

•  Microphysics:  Bnlk  schemes  ranging  from  simplihed  physics  snitable  for  mesoscale  modeling 

to  sophisticated  mixed-phase  physics  snitable  for  clond-resolving  modeling. 

•  Cumulus  parameterizations:  Adjustment  and  mass-flux  schemes  for  mesoscale  modeling  in- 

clnding  NWP. 

•  Surface  physics:  Mnlti-layer  land  surface  models  ranging  from  a  simple  thermal  model  to  fnll 

vegetation  and  soil  moistnre  models,  including  snow  cover  and  sea  ice. 

•  Planetary  boundary  layer  physics:  Turbulent  kinetic  energy  prediction  or  non-local  K  schemes. 

•  Atmospheric  radiation  physics:  Longwave  and  shortwave  schemes  with  mnltiple  spectral 

bands  and  a  simple  shortwave  scheme.  Clond  effects  and  surface  fluxes  are  included. 
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WRF-Var  System 


•  Incremental  formulation  of  the  model-space  cost  function. 

•  Quasi-Newton  or  conjugate  gradient  minimization  algorithms. 

•  Analysis  increments  on  un-staggered  Arakawa-A  grid. 

•  Representation  of  the  horizontal  component  of  background  error  B  via  recursive  hlters  (re¬ 

gional)  or  power  spectra  (global).  The  vertical  component  is  applied  through  projection 
onto  climatologically- averaged  eigenvectors  of  vertical  error.  Horizontal/vertical  errors  are 
non-separable  (horizontal  scales  vary  with  vertical  eigenvector). 

•  Background  cost  function  (Ji,)  preconditioning  via  a  control  variable  transform  U  dehned  as 

B  =  UU'^. 

•  Flexible  choice  of  background  error  model  and  control  variables. 

•  Climatological  background  error  covariances  estimated  via  either  the  NMC-method  of  aver¬ 

aged  forecast  differences  or  suitably  averaged  ensemble  perturbations. 

•  Unihed  3D-Var  (4D-Var  under  development),  global  and  regional,  multi-model  capability. 


WRF  Software  Framework 

•  Highly  modular,  single-source  code  for  maintainability. 

•  Portable  across  a  range  of  available  computing  platforms. 

•  Support  for  multiple  dynamics  solvers  and  physics  modules. 

•  Separation  of  scientihc  codes  from  parallelization  and  other  architecture-specihc  codes. 

•  Input/Output  Application  Program  Interface  (API)  enabling  various  external  packages  to  be 

installed  with  WRF,  hence  allowing  WRF  to  easily  support  various  data  formats. 

•  Efficient  execution  on  a  range  of  computing  platforms  (distributed  and  shared  memory,  vector 

and  scalar  types). 

•  Use  of  Earth  System  Modeling  Framework  (ESMF)  timing  package. 

•  Model  coupling  API  enabling  WRF  to  be  coupled  with  other  models  such  as  ocean,  and  land 

models. 


4 


Chapter  2 

Governing  Equations 


The  ARW  dynamics  solver  integrates  the  compressible,  nonhydrostatic  Euler  equations.  The 
equations  are  cast  in  flux  form  using  variables  that  have  conservation  properties,  following  the 
philosophy  of  Ooyama  (1990).  The  equations  are  formulated  using  a  terrain-following  mass 
vertical  coordinate  (Laprise,  1992).  In  this  chapter  we  define  the  vertical  coordinate  and  present 
the  flux  form  equations  in  Cartesian  space,  we  extend  the  equations  to  include  the  effects  of 
moisture  in  the  atmosphere,  and  we  further  augment  the  equations  to  include  projections  to  the 
sphere. 


2.1  Vertical  Coordinate  and  Variables 


The  ARW  equations  are  formulated  using  a 
terrain-following  hydrostatic-pressure  vertical  co¬ 
ordinate  denoted  by  rj  and  defined  as 

V  =  {Ph-Pht)/f^  where  ^^  =  Phs-Pht■  (2.1) 


0 


^ —  P^^  =  constant 


Ph  is  the  hydrostatic  component  of  the  pressure, 
and  phs  and  pht  refer  to  values  along  the  surface 
and  top  boundaries,  respectively.  The  coordinate 
definition  (2.1),  proposed  by  Laprise  (1992),  is 
the  traditional  a  coordinate  used  in  many  hydro¬ 
static  atmospheric  models,  rj  varies  from  a  value 
of  1  at  the  surface  to  0  at  the  upper  boundary  of 
the  model  domain  (Fig.  2.1).  This  vertical  coor¬ 
dinate  is  also  called  a  mass  vertical  coordinate. 

Since  p{x,  y)  represents  the  mass  per  unit  area 
within  the  column  in  the  model  domain  at  (x,  ?/), 
the  appropriate  flux  form  variables  are 

v  =  /iv  =  (f/,i/,iT),  n  =  pif],  e  =  pe.{2.2) 

V  =  (m,  V,  w)  are  the  covariant  velocities  in  the 
two  horizontal  and  vertical  directions,  respec¬ 
tively,  while  a;  =  ?7  is  the  contravariant  ‘vertical’ 


0.2 


Figure  2.1:  ARW  p  coordinate. 
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velocity.  6  is  the  potential  temperature.  Also  appearing  in  the  governing  equations  of  the  ARW 
are  the  non-conserved  variables  (j)  =  gz  (the  geopotential),  p  (pressure),  and  a  =  1/ p  (the  inverse 
density) . 


2.2  Flux-Form  Euler  Equations 

Using  the  variables  dehned  above,  the  flux-form  Euler  equations  can  be  written  as 


dtU  +  (V  ■  Yu)  -  a,(p0^)  +  a,(p0,) 

dtV  +  {Y  -Yv)-  dy{p(j)r,)  +  dy{p(py) 

dtW  +  (V  ■  Yw)  -  gidyP  -  p) 
dtQ  +  (V  ■  Ye) 

dtp+{Y-Y) 
dtct)  +  p-\Y-Y4>)-gW] 

along  with  the  diagnostic  relation  for  the  inverse  density 

dyCp  =  —ap, 

and  the  equation  of  state 

p  =  po{RdO/poa)'^. 


Fu 

(2.3) 

Fv 

(2.4) 

Fw 

(2.5) 

Fe 

(2.6) 

0 

(2.7) 

0 

(2.8) 

(2.9) 

(2.10) 

In  (2.3)  -  (2.10),  the  subscripts  x,  y  and  p  denote  differentiation. 


V  ■  Va  =  d^iUa)  +  dyiVa)  +  ^^(fla). 


and 

V  ■  Va  =  Ud^a  +  VdyU  +  Qd^a, 

where  a  represents  a  generic  variable.  7  =  Cp/cy  =  1.4  is  the  ratio  of  the  heat  capacities  for  dry 
air,  Rd  is  the  gas  constant  for  dry  air,  and  po  is  a  reference  pressure  (typically  10^  Pascals).  The 
right-hand-side  (RHS)  terms  Fu,  Fv,  Fw,  and  Fq  represent  forcing  terms  arising  from  model 
physics,  turbulent  mixing,  spherical  projections,  and  the  earth’s  rotation. 

The  prognostic  equations  (2.3)  ~  (2.8)  are  cast  in  conservative  form  except  for  (2.8)  which 
is  the  material  derivative  of  the  dehnition  of  the  geopotential.  (2.8)  could  be  cast  in  flux  form 
but  we  hnd  no  advantage  in  doing  so  since  pcf)  is  not  a  conserved  quantity.  We  could  also 
use  a  prognostic  pressure  equation  in  place  of  (2.8)  (see  Laprise,  1992),  but  pressure  is  not  a 
conserved  variable  and  we  could  not  use  a  pressure  equation  and  the  conservation  equation  for  0 
(2.6)  because  they  are  linearly  dependent.  Additionally,  prognostic  pressure  equations  have  the 
disadvantage  of  possessing  a  mass  divergence  term  multiplied  by  a  large  coefficient  (proportional 
to  the  sound  speed)  which  makes  spatial  and  temporal  discretization  problematic.  It  should  be 
noted  that  the  relation  for  the  hydrostatic  balance  (2.9)  does  not  represent  a  constraint  on  the 
solution,  rather  it  is  a  diagnostic  relation  that  formally  is  part  of  the  coordinate  dehnition.  In  the 
hydrostatic  counterpart  to  the  nonhydrostatic  equations,  (2.9)  replaces  the  vertical  momentum 
equation  (2.5)  and  it  becomes  a  constraint  on  the  solution. 
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2.3  Inclusion  of  Moisture 


In  formulating  the  moist  Euler  equations,  we  retain  the  coupling  of  dry  air  mass  to  the  prognostic 
variables,  and  we  retain  the  conservation  equation  for  dry  air  (2.7),  as  opposed  to  coupling  the 
variables  to  the  full  (moist)  air  mass  and  hence  introducing  source  terms  in  the  mass  conservation 
equation  (2.7).  Additionally,  we  define  the  coordinate  with  respect  to  the  dry-air  mass.  Based 
on  these  principles,  the  vertical  coordinate  can  be  written  as 

V  =  {Pdh  -  Pdht) / fJ^d  (2.11) 

where  jdd  represents  the  mass  of  the  dry  air  in  the  column  and  pdh  and  pdht  represent  the 
hydrostatic  pressure  of  the  dry  atmosphere  and  the  hydrostatic  pressure  at  the  top  of  the  dry 
atmosphere.  The  coupled  variables  are  defined  as 


V  =  /XrfV,  n  =  pdV,  0  =  Pdd. 

(2.12) 

With  these  dehnitions,  the  moist  Euler  equations  can  be  written  as 

dtU  +  (V  ■  Vu)r,  +  PdddxP  +  {a/ad)d^pdx4>  = 

Fu 

(2.13) 

dtV  +  (V  ■  Yv)rj  PdOidyP  +  {a/ad)dypdy4>  = 

Fv 

(2.14) 

dtW  +  (V  ■  Vw)y  -  g[{a/ad)dyp  -  pd]  = 

Fw 

(2.15) 

dtO  +  (V  ■  ve)y  = 

Fe 

(2.16) 

dtPd+{V-V)y  = 

0 

(2.17) 

dtcj>  +  Pd'[i^-V(j>)y-gW]  = 

0 

(2.18) 

dtQm  +  (V  ■  'Vqm)ri  = 

Fq^ 

(2.19) 

with  the  diagnostic  equation  for  dry  inverse  density 

drj(j)  OCdPd 

(2.20) 

and  the  diagnostic  relation  for  the  full  pressure  (vapor  plus  dry  air) 

P  =  Po{Rd9ra/  PoOldV 

(2.21) 

In  these  equations,  ad  is  the  inverse  density  of  the  dry  air  (1/pd)  and  a  is  the  inverse  density 
taking  into  account  the  full  parcel  density  a  =  ad{l  +  Qv  +  Qc  +  Qr  +  Qi  +  ■■■)~^  where  g*  are 
the  mixing  ratios  (mass  per  mass  of  dry  air)  for  water  vapor,  cloud,  rain,  ice,  etc.  Additionally, 

9m  ^(1  T  / Rd^Qv)  ~  ^(1  T  l-Slg^i),  and  Qm  PdQmi  Qm  Qvy  Qcj  Qii  ■ 

2.4  Map  Projections,  Coriolis  and  Curvature  Terms 

The  ARW  solver  currently  supports  three  projections  to  the  sphere —  the  Lambert  conformal, 
polar  stereographic,  and  Mercator  projections.  These  projections  are  described  in  Haltiner  and 
Williams  (1980).  These  projections,  and  the  ARW  implementation  of  the  map  factors,  assume 
that  the  map  factor  transformations  for  x  are  y  are  identical  at  a  given  point;  that  is,  the 
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transformation  is  isotropic.  Anisotropic  transformations,  snch  as  a  latitnde-longitnde  grid,  can 
be  accommodated  by  defining  separate  map  factors  for  the  x  and  y  transformations. 

In  the  ARW’s  compntational  space.  Ax  and  Ay  are  constants.  Orthogonal  projections  to 
the  sphere  reqnire  that  the  physical  distances  between  grid  points  in  the  projection  vary  with 
position  in  the  grid.  To  transform  the  governing  eqnations,  a  map  scale  factor  m  is  defined  as 
the  ratio  of  the  distance  in  compntational  space  to  the  corresponding  distance  on  the  earth’s 
snrface: 


(Ax,  Ay) 

m  =  — - ^ - -. 

distance  on  the  earth 

The  ARW  solver  inclndes  the  map-scale  factors  in  the  governing  eqnations  by  redefining  the 
moment  nm  variables  as 


(2.22) 


U  =  Hdu/m,  V  =  Hdv/m,  W  =  ydw/m,  hi  =  fidV/m. 


Using  these  redehned  momentnm  variables,  the  governing  eqnations,  inclnding  map  factors  and 


rotational  terms,  can  be  written  as 

dtU  +  m[da,{Uu)  +  dyiVu)]  +  dyiVtu)  +  ydOcd^P  +  {a/ad)dypda:(j)  =  Fu  (2.23) 

dtV  +  m[dj,{Uv)  +  dy{Vv)]  +  dy{Qv)  pdOidyP  {a/ad)drjpdy4>  =  Fy  (2.24) 

dtW  m[d^{Uw)  dy{Vw)]  dy{^lw)  -  m~^g[{a/ad)dyp  -  pd]  =  Fw  (2.25) 

dtQ  +  m^[d^{ue)  +  dy{ve)]  +  mdy{ne)  =  Fe  (2.26) 

dtPd  +  m‘^[Ux  +  Vy]+mdri{n)  =  0  (2.27) 

dtcj)  +  +  V4>y)  +  rnVL(l)y  —  gW]  =  0  (2.28) 

dtQm  +  rn^[dx{Uq^)  +  dyiV  q^)]  +  mdyiVLqyn)  =  Fq^,  (2.29) 

and,  for  completeness,  the  diagnostic  relation  for  the  dry  inverse  density 

dy(j)  =  -adPd,  (2.30) 

and  the  diagnostic  eqnation  for  fnll  pressnre  (vapor  pins  dry  air) 

p  =  PoiRdOm/p^adV .  (2.31) 


The  right-hand-side  terms  of  the  momentnm  eqnations  (2.23)  -  (2.25)  contain  the  Coriolis 
and  curvatnre  terms  along  with  mixing  terms  and  physical  forcings.  The  Coriolis  and  cnrvatnre 
terms  can  be  written  as  follows: 


^  ^  dm  dm\^^  uW 

Fu„.^  ^+{f  +  u—-v^jV-eWco.ar-  — 

^  „  dm  dm\^^  vW 

n„,  =  -I  /  +  jc'  +  +  — 


Fw^or  =  +e(U  cosa^  -  V  sina^) 


uU  +  vV 


(2.32) 

(2.33) 

(2.34) 


where  is  the  local  rotation  angle  between  the  y-axis  and  the  meridians,  ^|J  is  the  latitnde, 
/  =  2r2eSin'0,  e  =  2VteCos'ip,  is  the  angnlar  rotation  rate  of  the  earth,  and  Xg  is  the  radins  of 


the  earth.  In  this  formulation  we  have  approximated  the  radial  distance  from  the  center  of  the 
earth  as  the  mean  earth  radius  Ve,  and  we  have  not  taken  into  account  the  change  in  horizontal 
grid  distance  as  a  function  of  the  radius.  The  terms  containing  m  are  the  horizontal  curvature 
terms,  those  containing  relate  to  vertical  (earth-surface)  curvature,  and  those  with  e  and  / 
are  the  Coriolis  force.  In  idealized  cases,  the  map  scale  factor  m  =  1,  /  is  often  taken  to  be 
constant,  and  e  =  0. 

2.5  Perturbation  Form  of  the  Governing  Equations 

Before  constructing  the  discrete  solver,  it  is  advantageous  to  recast  the  governing  equations 
using  perturbation  variables  so  as  to  reduce  truncation  errors  in  the  horizontal  pressure  gradient 
calculations  in  the  discrete  solver,  in  addition  to  reducing  machine  rounding  errors  in  the  vertical 
pressure  gradient  and  buoyancy  calculations.  For  this  purpose,  new  variables  are  dehned  as 
perturbations  from  a  hydrostatically-balanced  reference  state,  and  we  dehne  reference  state 
variables  (denoted  by  overbars)  that  are  a  function  of  height  only  and  that  satisfy  the  governing 
equations  for  an  atmosphere  at  rest.  That  is,  the  reference  state  is  in  hydrostatic  balance  and 
is  strictly  only  a  function  of  z.  In  this  manner,  p  =  p{z)  +  p' ,  (j)  =  0(z)  +  (p',  a  =  a{z)  +  a', 
and  pd  =  Jid{x,y)  +  p'^.  Because  the  p  coordinate  surfaces  are  generally  not  horizontal,  the 
reference  prohles  p,  0,  and  a  are  functions  of  {x,y,p).  The  hydrostatically  balanced  portion 
of  the  pressure  gradients  in  the  reference  sounding  can  be  removed  without  approximation  to 
the  equations  using  these  perturbation  variables.  The  momentum  equations  (2.23)  -  (2.25)  are 
written  as 

dtU  m[dx{Uu)  +  dy{Vu)\  +  dy{Vtu)  +  {pdad^p  +  PdCt'd^p) 

+  {a/ad){pddx(j)'  +  dyp'd^cj)  -  PdO^cf))  =  Fu  (2.35) 

dtV  +  m[d^{Uv)  +  dy{Vv)]  +  drj{D.v)  +  {pdadyp  +  pdadyp) 

+  {a/ad){pddy(j)'  +  dypdycj)  -  p'ddycf))  =  Fy  (2.36) 

dfW  +  m[dx{Uw)  +  dy{Vw)]  +  dy{Qw) 

-m~^g{a/ ad)  [dyp  -  pd{qv  +  qc  +  qr)]  +  rn~^p^g  =  Fw,  (2.37) 

and  the  mass  conservation  equation  (2.27)  and  geopotential  equation  (2.28)  become 

dtp'd  +  m^[dxU  +  dyV]  -|-  mdyVt  =  0  (2.38) 

dtcf)'  +  p'2^[w?{U(f)x  +  y  (f>y)  +  rnilcpy  —  gW]  =  0.  (2.39) 

Remaining  unchanged  are  the  conservation  equations  for  the  potential  temperature  and  scalars 

ds  +  m^[dx{ue)  +  dyive)]  +  mdyiyte)  =  f©  (2.40) 

dtQni  +  rn^[dx{Uqyn)  +  dyiV  qm)]  +  mdyiVLqm)  =  Fq^.  (2.41) 

In  the  perturbation  system  the  hydrostatic  relation  (2.30)  becomes 

^?70  ^d^^d'  (2.42) 

Equations  (2.35)  -  (2.41),  together  with  the  equation  of  state  (2.21),  represent  the  equations 
solved  in  the  ARW.  The  RHS  terms  in  these  equations  include  the  Coriolis  terms  (2.32)  - 
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(2.34),  mixing  terms  (described  in  Chapter  4),  and  parameterized  physics  (described  in  Chapter 
8).  Also  note  that  the  equation  of  state  (2.21)  cannot  be  written  in  perturbation  form  because 
of  the  exponent  in  the  expression.  For  small  perturbation  simulations,  accuracy  for  perturbation 
variables  can  be  maintained  by  linearizing  (2.21)  for  the  perturbation  variables. 
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Chapter  3 

Model  Discretization 


3.1  Temporal  Discretization 

The  ARW  solver  uses  a  time-split  integration  scheme.  Generally  speaking,  slow  or  low-frequency 
(meteorologically  significant)  modes  are  integrated  using  a  third-order  Runge-Kutta  (RK3)  time 
integration  scheme,  while  the  high-frequency  acoustic  modes  are  integrated  over  smaller  time 
steps  to  maintain  numerical  stability.  The  horizontally  propagating  acoustic  modes  (including 
the  external  mode  present  in  the  mass-coordinate  equations  using  a  constant-pressure  upper 
boundary  condition)  are  integrated  using  a  forward-backward  time  integration  scheme,  and  ver¬ 
tically  propagating  acoustic  modes  and  buoyancy  oscillations  are  integrated  using  a  vertically 
implicit  scheme  (using  the  acoustic  time  step).  The  time-split  integration  is  similar  to  that  hrst 
developed  by  Klemp  and  Wilhelmson  (1978)  and  analyzed  by  Skamarock  and  Klemp  (1992). 
The  time-split  RK3  scheme  is  described  in  general  terms  in  Wicker  and  Skamarock  (2002). 
The  primary  differences  between  the  descriptions  found  in  the  references  and  the  ARW  imple¬ 
mentation  are  associated  with  our  use  of  the  mass  vertical  coordinate  and  a  flux-form  set  of 
equations,  along  with  our  use  of  perturbation  variables  for  the  acoustic  component  of  the  time- 
split  integration.  The  acoustic-mode  integration  is  cast  in  the  form  of  a  correction  to  the  RK3 
integration. 

3.1.1  Runge-Kutta  Time  Integration  Scheme 

The  RK3  scheme,  described  in  Wicker  and  Skamarock  (2002),  integrates  a  set  of  ordinary 
differential  equations  using  a  predictor-corrector  formulation.  Defining  the  prognostic  variables 
in  the  ARW  solver  as  $  =  {U,V,W,Q,(f)' ,  fi' ,Qm)  and  the  model  equations  as  =  !?($),  the 
RK3  integration  takes  the  form  of  3  steps  to  advance  a  solution  ^(t)  to  +  At): 

$*  =  $*  + —R(d)‘)  (3.1) 

(3.) 

^t+At  =  ^  AtR{^**)  (3.3) 

where  At  is  the  time  step  for  the  low-frequency  modes  {the  model  time  step).  In  (3.1)  -  (3.3), 
superscripts  denote  time  levels.  This  scheme  is  not  a  true  Runge-Kutta  scheme  per  se  because. 
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while  it  is  third-order  accurate  for  linear  equations,  it  is  only  second-order  accurate  for  nonlinear 
equations.  With  respect  to  the  ARW  equations,  the  time  derivatives  are  the  partial  time 
derivatives  (the  leftmost  terms)  in  equations  (2.35)  -  (2.41),  and  i?(*h)  are  the  remaining  terms 
in  (2.35)  -  (2.41). 


3.1.2  Acoustic  Integration 

The  high-frequency  but  meteorologically  insignihcant  acoustic  modes  would  severely  limit  the 
RK3  time  step  At  in  (3.1)  -  (3.3).  To  circumvent  this  time  step  limitation  we  use  the  approach 
described  in  Wicker  and  Skamarock  (2002).  Additionally,  to  increase  the  accuracy  of  the  split¬ 
ting,  we  integrate  a  perturbation  form  of  the  governing  equations  using  smaller  acoustic  time 
steps  within  the  RK3  large-time-step  sequence.  To  form  the  perturbation  equations  for  the  RK3 
time-split  acoustic  integration,  we  dehne  small  time  step  variables  that  are  deviations  from  the 
most  recent  RK3  predictor  (denoted  by  the  superscript  t*  and  representing  either  <!>*,  $*,  or  4)** 
in  (3.1)  -  (3.3)): 

v"  =  v-v**,  =  0"  =  0-0^\ 

in  1/  i/t*  n  II  t*  //  /  it* 

0=0-0  ,  a^  =  a^-a^  ,  =  /i^  -  /r  ^  . 

The  hydrostatic  relation  (i.e.,  the  vertical  coordinate  dehnition)  becomes 


Additionally,  we  also  introduce  a  version  of  the  equation  of  state  that  is  linearized  about  f*. 


(3.5) 


where  cf  =  7P**a^*  is  the  square  of  the  sound  speed.  The  linearized  state  equation  (3.5)  and 
the  vertical  coordinate  dehnition  (3.4)  are  used  to  cast  the  vertical  pressure  gradient  in  (2.37) 
in  terms  of  the  model’s  prognostic  variables.  By  combining  (3.5)  and  (3.4),  the  vertical  pressure 
gradient  can  be  expressed  as 


5,p"  =  (3.6) 

where  C  =  This  linearization  about  the  most  recent  large  time  step  should  be  highly 

accurate  over  the  time  interval  of  the  several  small  time  steps. 

These  variables  along  with  (3.6)  are  substituted  into  the  prognostic  equations  (2.35)  -  (2.41) 
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and  lead  to  the  acoustic  time-step  equations: 


5rU''  +  {i/*d^p)a"'^  +  {a/ad)[f/* -f  {d^(t/*){dnp'  -  p'Y] 

6rV"  +  Y*a^*dyp”^  +  {Y*dyp)a"^  +  {a/ad)[Y*dyY'^  +  {dyYYidrjP"  -  pY] 

5rYd  +  mYd^U”  +  dyV'Y^^"  + 

5,0"  +  mYdYu”o^*)  +  dy{y''e^YY^^"  +  mdYQ!'^^^^e^Y 


YW  -  m-^g 


{a/ady*dYCdyY')  +  9y 


Y  Q" 


Yd 


1 


5,0"  +  -  gW'Y 


R/ 

Re‘‘ 

Rw‘' 

R/. 


(3.7) 

(3.8) 

(3.9) 

(3.10) 

(3.11) 

(3.12) 


The  RHS  terms  in  (3.7)  -  (3.12)  are  hxed  for  the  acoustic  steps  that  comprise  the  time  integration 
of  each  RK3  sub-step  (i.e.,  (3.1)  -  (3.3)),  and  are  given  by 


RY  =  -  m[dYUu)  +  dy{Vu)]  dy{D.u)  -  {pdCtd^p  -  pdYd^p) 


-  {a/ad){pddxY  -  dypdxcf)  +  Yd^xY)  +  Fu  (3.13) 

HY  =  -  m[dx{Uv)  +  dy{Vv)]  +  dy{^lv)  -  {pdadyp'  -  pdYdyp) 

-  {a/ad){pddyY  -  dypdycf)  Yd^yY)  +  Fy  (3-14) 

=  —  nY[dxU  +  dyV]  -|-  mdyVt  (3.15) 

rY  =  -  mYdxiue)  +  dyive)]  -  mdyiyte)  +  f©  (3.16) 

RY/  =  —  m[dx{Uw)  +  dyiVw)]  +  dy{Q,w) 

-h  m~^g{a/ad)[dyp  fldYv  +  qc  +  qY]  -  Yd9  +  Fw  (3.17) 

^5)*  =  -  pY['>^‘^iU(l>x  +  V(j)y) +m^l(j)y- gW],  (3.18) 


where  all  variables  in  (3.13)  ~  (3.18)  are  evaluated  at  time  t*  (i.e.,  using  <h*,  $*,  or  $**  for  the 
appropriate  RK3  sub-step  in  (3.1)  -  (3.3)).  Equations  (3.7)  -  (3.12)  utilize  the  discrete  acoustic 
time-step  operator 


r+Ar 


5,a  = 


At 


where  At  is  the  acoustic  time  step,  and  an  acoustic  time-step  averaging  operator 


\  (3 


Y  =  + 


1-/3 


(3.19) 


where  /3  is  a  user-specihed  parameter  (see  Section  4.2.3). 

The  integration  over  the  acoustic  time  steps  proceeds  as  follows.  Beginning  with  the  small 
time-step  variables  at  time  r,  (3.7)  and  (3.8)  are  stepped  forward  to  obtain  and 

Both  and  are  then  calculated  from  (3.9).  This  is  accomplished  by  Erst  integrating 

(3.9)  vertically  from  the  surface  to  the  material  surface  at  the  top  of  the  domain,  which  removes 
the  dyVt''  term  such  that 


5rPd  =  rY  /  [dxU"  +  dyV'Y^^^dp. 


(3.20) 
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After  computing  from  (3.20),  is  recovered  by  using  (3.9)  to  integrate  the  dr^Vt" 

term  vertically  using  the  lower  boundary  condition  =  0  at  the  surface.  Equation  (3.10) 
is  then  stepped  forward  to  calculate  0"'^+^'^.  Equations  (3.11)  and  (3.12)  are  combined  to 
form  a  vertically  implicit  equation  that  is  solved  for  subject  to  the  boundary  condition 

W"  =  V"  ■  Vh  at  the  surface  {z  =  h{x,y))  and  p'  =  0  along  the  model  top.  is  then 

obtained  from  (3.12),  and  and  are  recovered  from  (3.5)  and  (3.4). 

3.1.3  Full  Time-Split  Integration  Sequence 

The  time-split  RK3  integration  technique  is  summarized  below.  It  consists  of  two  primary 
loops —  an  outer  loop  for  the  large-time-step  Runge-Kutta  integration,  and  an  inner  loop  for 
the  acoustic  mode  integration. 


Begin  Time  Step 

Begin  RK3  Loop:  Steps  1,  2,  and  3 

(1)  If  RK3  step  1,  compute  and  store  Fq> 

(i.e.,  physics  tendencies  for  RK3  step,  including  mixing). 

(2)  Compute  R%^  (3.13)-(3.18) 

Begin  Acoustic  Step  Loop:  Steps  1  — n, 

RK3  step  1,  n  =  1,  At  =  At/3; 

RK3  step  2,  n  =  71^/2,  At  =  At/ug; 

RK3  step  3,  n  =  Ug,  At  =  At/ug. 

(3)  Advance  horizontal  momentum,  (3.7)  and  (3.8) 

(4)  Advance  pd  (3.9)  and  compute  then  advance  0  (3.10) 

(5)  Advance  W  and  0  (3.11)  and  (3.12) 

(6)  Diagnose  p"  and  a"  using  (3.5)  and  (3.4) 

End  Acoustic  Step  Loop 

(7)  Scalar  transport:  Advance  scalars  (2.41) 

over  RK3  substep  (3.1),  (3.2)  or  (3.3) 

(using  mass  fluxes  U,  V  and  D  time-averaged  over  the  acoustic  steps). 

(8)  Compute  p'  and  a'  using  updated  prognostic  variables  in  (2.31)  and  (2.42) 

End  RK3  Loop 

(9)  Compute  non-RK3  physics  (currently  microphysics),  advance  variables. 

End  Time  Step 


Figure  3.1:  Time  step  integration  sequence.  Here  n  represents  the  number  of  acoustic  time 
steps  for  a  given  substep  of  the  RK3  integration,  and  Ug  is  the  ratio  of  the  RK3  time  step  to 
the  acoustic  time  step  for  the  second  and  third  RK3  substeps. 

In  the  RK3  scheme,  physics  can  be  integrated  within  the  RK3  time  integration  (using  a  time 
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forward  step,  i.e.,  step  (1)  in  Fig.  3.1,  or  the  RK3  time  integration  if  higher  temporal  accuracy 
is  desired,  i.e.,  in  step  (2) —  implying  a  physics  evaluation  every  RK3  substep)  or  external  to  it 
using  additive  timesplitting,  i.e.,  step  (9). 

Within  the  acoustic  integration,  the  acoustic  time  step  At  is  specihed  by  the  user  through 
the  choice  of  Ug  (see  Section  3.3.2).  Within  the  hrst  RK3  substep,  however,  a  single  acoustic 
time  step  is  used  to  advance  the  solution  regardless  of  Ug.  Within  the  full  RK3- acoustic  timesplit 
integration,  this  modihed  acoustic  time  step  does  not  impose  any  additional  stability  constraints 
(see  Wicker  and  Skamarock,  2002). 

The  major  costs  in  the  model  arise  from  the  evaluation  of  the  right  hand  side  terms  i?**  in 
(3.7)  -  (3.12).  The  efficiency  of  the  RK3  timesplit  scheme  arises  from  the  fact  that  the  RK3 
time  step  At  is  much  larger  than  the  acoustic  time  step  At,  hence  the  most  costly  evaluations 
are  only  performed  in  the  less-frequent  RK3  steps. 


3.1.4  Diabatic  Forcing 

Within  the  RK3  integration  sequence  outlined  in  Fig.  3.1,  the  RHS  term  Rq  in  the  thermo¬ 
dynamic  equation  (3.10)  contains  contributions  from  the  diabatic  physics  tendencies  that  are 
computed  in  step  (1)  at  the  beginning  of  the  hrst  RK3  step.  This  diabatic  forcing  is  integrated 
within  the  acoustic  steps  (specihcally,  in  step  4  in  the  time  integration  sequence  shown  in  Fig. 
3.1).  Additional  diabatic  contributions  are  integrated  in  an  additive-time-split  manner  in  step 
(9)  after  the  RK3  update  is  complete.  Thus,  the  diabatic  forcing  computed  in  step  (9)  (the 
microphysics  in  the  current  release  of  the  ARW)  does  not  appear  in  Rq  from  (3.10)  in  the 
acoustic  integration.  We  have  discovered  that  this  time  splitting  can  excite  acoustic  waves  and 
can  give  rise  to  noise  in  the  solutions  for  some  applications.  Note  that  the  non-RK3  physics  are 
integrated  in  step  (9)  because  balances  produced  in  the  physics  are  required  at  the  end  of  the 
time  step  (e.g.,  the  saturation  adjustment  in  the  microphysics).  So  while  moving  these  non-RK3 
physics  into  step  (1)  would  eliminate  the  noise,  the  balances  produced  by  these  physics  would 
be  altered. 

We  have  found  that  the  excitation  of  the  acoustic  modes  can  be  circumvented  while  leaving 
the  non-RK3  physics  in  step  (9)  by  using  the  following  procedure  that  is  implemented  in  the 
ARW.  In  step  (1)  of  the  integration  procedure  (Fig.  3.1),  an  estimate  of  the  diabatic  forcing 
of  0  arising  from  the  non-RK3  physics  in  step  (9)  is  included  in  the  diabatic  forcing  term 
Rq  in  (3.10)  (which  is  advanced  in  step  4).  This  estimated  diabatic  forcing  is  then  removed 
from  the  updated  0  after  the  RK3  integration  is  complete  and  before  the  evaluation  of  the 
non-RK3  physics  in  step  (9).  We  use  the  diabatic  forcing  from  the  previous  time  step  as  the 
estimated  forcing;  hence  this  procedure  results  in  few  additional  computations  outside  of  saving 
the  diabatic  forcing  between  time  steps. 


3.1.5  Hydrostatic  Option 

A  hydrostatic  option  is  available  in  the  ARW  solver.  The  time-split  RK3  integration  technique 
summarized  in  Fig.  3.1  is  retained,  including  the  acoustic  step  loop.  Steps  (5)  and  (6)  in  the 
acoustic-step  loop,  where  W  and  cf)  are  advanced  and  p"  and  a"  are  diagnosed,  are  replaced  by 
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(1),  the  diagnosis  of  the  hydrostatic  pressure  using  the  dehnition  of  the  vertical  coordinate 

^rjPh  ^  ^  Qm)  Pdi 

followed  by  (2),  the  diagnosis  of  ad  using  the  equation  of  state  (2.31)  and  the  prognosed  9,  and 
(3),  the  diagnosis  of  the  geopotential  using  the  hydrostatic  equation 

Pd^d  Pd^d- 

The  vertical  velocity  w  can  be  diagnosed  from  the  geopotential  equation,  but  it  is  not  needed  in 
the  solution  procedure.  The  acoustic  step  loop  advances  gravity  waves,  including  the  external 
mode,  when  the  hydrostatic  option  is  used;  there  are  no  horizontally  propagating  acoustic  modes 
in  this  hydrostatic  system. 


3.2  Spatial  Discretization 

The  spatial  discretization  in  the  ARW  solver  uses  a  C  grid  staggering  for  the  variables  as  shown 
in  Fig.  3.2.  That  is,  normal  velocities  are  staggered  one-half  grid  length  from  the  thermodynamic 
variables.  The  variable  indices,  {i,j)  for  the  horizontal  plane  and  {i,k)  for  the  vertical  plane, 
indicate  variable  locations  where  {x,y,ri)  =  {iAx,  jAy,  kArj).  We  will  denote  the  points  where 
9  is  located  as  being  mass  points,  and  likewise  we  will  denote  locations  where  u,  v,  and  w  are 
dehned  as  u  points,  v  points,  and  w  points,  respectively.  Not  shown  in  Fig.  3.2  are  the  column 
mass  /X,  dehned  at  the  {i,j)  points  (mass  points)  on  the  discrete  grid,  the  geopotential  0  that 
is  dehned  at  the  w  points,  and  the  moisture  variables  qm  are  dehned  at  the  mass  points.  The 
diagnostic  variables  used  in  the  model,  the  pressure  p  and  inverse  density  a,  are  computed  at 
mass  points.  The  grid  lengths  Ax  and  Ay  are  constants  in  the  model  formulation;  changes  in 
the  physical  grid  lengths  associated  with  the  various  projections  to  the  sphere  are  accounted 
for  using  the  map  factors  introduced  in  Section  2.4.  The  vertical  grid  length  Ap  is  not  a  hxed 
constant;  it  is  specihed  in  the  initialization.  The  user  is  free  to  specify  the  p  values  of  the  model 
levels  subject  to  the  constraint  that  ?7  =  1  at  the  surface,  ?7  =  0  at  the  model  top,  and  p  decreases 
monotonically  between  the  surface  and  model  top.  Using  these  grid  and  variable  dehnitions,  we 
can  dehne  the  spatial  discretization  for  the  ARW  solver. 

3.2.1  Acoustic  Step  Equations 

We  begin  by  dehning  the  column-mass-coupled  variables  relative  to  the  uncoupled  variables. 
The  vertical  velocity  is  staggered  only  in  fc,  so  it  can  be  coupled  directly  to  the  column  mass 
with  no  averaging  or  interpolation.  The  horizontal  velocities  are  horizontally  staggered  relative 
to  the  column  mass  such  that  the  continuous  variables  are  represented  discretely  as 

JJ  _  PdU  ^  JI^U  _  PdV  ^  Jj^v 

m  m  m 

where  the  discrete  operator  denotes  a  linear  interpolation  operator.  The  grid  lengths  Ax  and 
Ay  are  constant,  hence  in  this  case  the  operator  reduces  to  =  (ai+i/2  -l-  ai-ii2)/2. 


16 


ri 


Ay 


horizontal  grid  vertical  grid 

Figure  3.2:  Horizontal  and  vertical  grids  of  the  ARW 

Using  these  dehnitions,  we  can  write  the  spatially  discrete  acoustic  step  equations 
(3.7)  -  (3.12)  as 


drU”  +  jJ*  a**  6xp"^  +  (/i**  ^xP)ci'‘ 


+  {6x4>^*\5^p''^  '  -  p"y]  =  Rl 


5rV"  +  y^a**^6yp"^  +  6yp)a”'’ 


+y/y^[W^6yy^  +  (5,r  )( V"  -  J^Y]  =  Rv 

5rPd  +  mYSxU"  +  6yV''Y+^^  + 

drQ"  +  m^yU'W")  +  6y{V"¥^)Y+^^  +  m6y{n"^+^W)  =  Rq' 


5rW''  -  m-^g 


[a/ ay  Yicdy)  +  Y 


cj  0" 


=R 


w 


dy  +  -  gW'Y  =  R/*, 

Yd 


where  the  discrete  operator 


6xCl  —  Ax  {cii+l/2  ~  o.i-1/2) 


(3.21) 

(3.22) 

(3.23) 

(3.24) 

(3.25) 

(3.26) 


(3.27) 


with  the  operators  5y  and  5y  similarily  dehned.  Additionally,  the  operator  aP  is  a  vertical 
interpolation  operator.  Using  the  notation  given  for  the  vertically  stretched  grid  depicted  in 
Fig.  3.2,  it  is  dehned  as 


(R\k+l/2 


^  f 


Apfc+i  \ 

Apfc+1/2  v' 


(3.28) 
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The  operator  vertically  interpolates  variables  on  mass  levels  k  to  the  w  levels  (fc  +  |).  It  should 
be  noted  that  the  vertical  grid  is  dehned  such  that  vertical  interpolation  from  w  levels  to  mass 
levels  reduces  to  =  {ak+1/2  +  afc-i/2)/2  (see  Fig.  3.2). 

The  RHS  terms  in  the  discrete  acoustic  step  equations  for  momentum  (3.21),  (3.22)  and 
(3.25)  are  discretized  as 


Ru 


Rv 


+  +  advection  +  mixing  +  physics 

-  {jM^a^6yP  -  -  {a/adf  {jld^6y(j)''^  -  5yp'^  ^ dycjP  +  5y(tP) 


+  +  advection  +  mixing  +  physics 

=  m-^g^a/adflSyp'  +  pdPFP]  - 

+  FVcor  +  advection  +  mixing  +  buoyancy  +  physics. 


(3.29) 

(3.30) 

(3.31) 


3.2.2  Coriolis  and  Curvature  Terms 


The  terms  Fu^^^,  Fy^^^,  and  FVcor  (3.29)  -  (3.31)  represent  Coriolis  and  curvature  effects 
in  the  equations.  These  terms  in  continuous  form  are  given  in  (2.32)  -  (2.34).  Their  spatial 
discretization  is 

=  +{T  +  u^5ym  -  cbso^^  - 

^  e 

Fy_  =  -{T  +  u^6ym-vy6,m^)  T" "  +  e  W""  ^ 

Fw^or  =  COS  ar  -  sin  a r)  + 


-r^rxy 


Here  the  operators  ()''  =  ()  ,  and  likewise  for  ()  ^  and  () 


3.2.3  Advection 

The  advection  terms  in  the  ARW  solver  are  in  the  form  of  a  flux  divergence  and  are  a  subset  of 


the  RHS  terms  in  equations  (3.13)  -  (3.18): 

=  -  m[da,{Uu)  +  dy{Vu)]  +  dy{nu)  (3.32) 

Rvadv  ^  ~  rn[dx{Uv)  +  dy{Vv)]  +  dy{Q,v)  (3.33) 

Rtad.  =  -  +  ^y]  +  (3-34) 

=  -  m‘^[d^{Ue)  +  dy{Ve)]  -  mdy{ne)  (3.35) 

Rpdv  ^  ~  m[dx{Uw)  +  dy{Vw)]  +  dy{Qw)  (3.36) 

R^ad.  =  -  +  V(j)y)+  mn(j)y].  (3.37) 

For  the  mass  conservation  equation,  the  flux  divergence  is  discretized  using  a  2nd-order  centered 
approximation: 

=  -m^[6^U  +  6yVf  +  rndykl^*.  (3.38) 
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In  the  current  version  of  the  ARW,  the  advection  of  vector  quantities  (momentum)  and  scalars 
is  performed  using  the  RK3  time  integration  as  outlined  in  Fig.  3.1.  The  spatial  discretization 
used  in  this  approach  is  outlined  in  the  next  section.  For  many  applications  it  is  desirable  to 
use  positive  dehnite  or  monotonic  advection  schemes  for  scalar  transport.  In  the  next  major 
release  of  the  ARW  we  will  be  including  a  forward-in-time  scheme  for  scalar  transport  that  has 
positive  dehnite  and  monotonic  options.  We  describe  that  scheme  in  the  section  following  the 
description  of  the  RK3  advection. 

RK3  Advection 

2nd  |;P2;ough  6*^  order  accurate  spatial  discretizations  of  the  hux  divergence  are  available  in  the 
ARW  for  momentum,  scalars  and  geopotential  using  the  RK3  time-integration  scheme  (scalar 
advection  option  1,  step  7  in  the  time-split  integration  sequence  in  Fig.  3.1).  The  discrete 
operators  can  be  illustrated  by  considering  the  hux  divergence  equation  for  a  scalar  q  in  its 
discrete  form: 

(3.39) 

As  in  the  pressure  gradient  discretization,  the  discrete  operator  is  dehned  as 

=  Ax-'  [(i7r“‘''')*+i/2  -  (t^r“'^'')i-i/2] .  (3.40) 

The  diherent  order  advection  schemes  correspond  to  diherent  dehnitions  for  the  operator  _ 
The  even  order  operators  (2"''^,  4'^,  and  are 

2”“^  order:  (g^“'‘")i_i/2  =  +  ?i-i) 

7 

d''*  order:  =  —(?*  +  Qi-i) 
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order:  (g^“'*")i_i/2  =  +  ?i-i) 

dU 

and  the  odd  order  operators  (3^'^  and  5'^)  are 

3^-^  order:  (r“''’^),_i/2  =  (r“'^’')fi/2 

+  sign(f/)^  [(gi+i  -  qi_2)  -  3{qi  -  g^.i)] 

5^^  order:  (g^“''”)i-i/2  =  (g''“''’')tv2 

-  sign(t/)^  [(gi+2  -  qi-s)  -  5(gi+i  -  qi-2)  +  10(gi  -  g^-i)] . 

The  even-order  advection  operators  are  spatially  centered  and  thus  contain  no  implicit  dif¬ 
fusion  outside  of  the  dihusion  inherent  in  the  RK3  time  integration.  The  odd-order  schemes  are 
upwind-biased,  and  the  spatial  discretization  is  inherently  dihusive.  As  is  evident  in  their  for¬ 
mulation,  the  odd-order  schemes  are  comprised  of  the  next  higher  (even)  order  centered  scheme 
plus  an  upwind  term  that,  for  a  constant  transport  mass  hux,  is  a  dihusion  term  of  that  next 
higher  (even)  order  with  a  hyper- viscosity  proportional  to  the  Courant  number  {Cr).  Further 
details  concerning  RK3  advection  can  be  found  in  Wicker  and  Skamarock  (2002) 


“  +  ^*-2) 

2  /  N  1  / 

—  —{qi+1  +  (li-2)  +  '^{Qi+2  +  Qi-3)^ 
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Forward-In-Time  Scalar  Advection 


A  forward-in-time  scalar  advection  scheme  having  monotonic  and  positive  dehnite  options  will 
be  available  in  the  next  major  release  of  the  ARW.  This  option  entails  bypassing  step  (7)  in  the 
time-split  integration  sequence  in  Fig.  3.1,  and  adding  a  single  advection  evaluation  before  step 
(9),  after  the  end  of  the  RK3  loop  in  the  integration  sequence.  The  new  advection  algorithm  is 
patterned  after  Easter  (1993)  and  can  be  written  as  follows 


(Mdq)*  =  il^dqY  +  m^FxiqY 

(3.41) 

^^d  =  ^^d  +  ^^^-'^FxiI) 

(3.42) 

Q*  =  (mT/M* 

(3.43) 

if^dq)**  =  ifJ^dq)*  +m'^Fy{q*) 

(3.44) 

(3.45) 

q**  =  M**/M** 

(3.46) 

{fidqY^^^  =  {fidqr+mFy{qn 

(3.47) 

(3.48) 

(3.49) 

(3.49)  the  operator 

Fx{q)  =  -AtAx~^{{U^q)i+i/2  -  {U^q)i-i/2), 

(3.50) 

with  similar  dehnitions  for  Fy  and  and  /  is  a  vector  with  all  values  equal  to  1.  The  discrete 
mass  continuity  equation  ((3.42)-t-(3.45)-l-(3.48))  can  be  written  as 

.  (3.51) 


It  can  easily  be  seen  that  the  scheme  (3.41)  -  (3.49)  collapses  to  (3.51)  for  q  =  I,  and  hence  it  is 
consistent.  The  mass  fluxes  U  ,  V  ,  and  ,  represent  time- averaged  values  where  the  averaging 
is  performed  over  the  hnal  RK3  small-time-step  cycle.  Hence,  the  mass  conservation  equation 
(3.51)  produces  the  mass  conservation  equation  integrated  within  the  RK3  scheme.  Equation 
(3.51)  is  re-integrated  within  the  time-split  transport  scheme  only  because  a  consistent  column 
mass  prf  is  needed  on  the  sub-steps  to  retrieve  q  from  the  prognostic  variable 

Any  forward-in-time  flux  operator  that  is  stable  for  Cr  <  1  can  be  used  to  evaluate  the 
operators  F^^q),  Fy{q),  or  Fy{q).  A  scheme  based  on  the  Piecewise  Parabolic  Method  (PPM) 
has  been  implemented  in  the  ARW.  It  is  the  non-monotonized  PPM  advection  described  in 
Carpenter  et  al.  (1990),  and  it  provides  the  fluxes  {U  q)i±i  used  in  (3.50)  in  the  flux  divergence 
operators.  Dehning  g*  as  the  control-volume  average  mixing  ratio  for  cell  i,  zone  edge  values 
qadv  dehned  as 

=  (7('?j+i  +  Qi)  ~  {Qi+2  +  5'i-i))/12.  (3.52) 

The  flux  through  the  i  +  \  face  can  be  written  as 


i-\-  ■ 


[gf^l  -  C'r(g“fl  -  g,)  -  C'r(l  -  C'r)(g“fl  -  2g,  +  g“fl)] 


(3.53) 
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for  Cr  >  0,  and 

(u‘qU^  =  +  Cr(l  +  Cr)(,^”.  -  2,.+.  +  qffl )]  (3.54) 

for  Cr  <  0,  where  Cr  =  u^j_iAt/ Ax.  In  addition,  to  maintain  second-order  accuracy  for  the 
time-split  forward-in-time  advection  scheme,  the  sequence  (3.41)  -  (3.49)  is  reversed  every  other 
time  step;  that  is,  the  flux  divergence  in  rj  is  computed  hrst  followed  by  y  and  then  x.  Finally, 
the  scheme  has  also  been  augmented  for  Courant  numbers  greater  than  one  using  the  ID  “semi- 
Lagrangian”  flux  calculation  described  in  Lin  and  Rood  (1996).  More  information  about  this 
scheme  and  a  description  of  the  limiters  can  be  found  in  Skamarock  (2005). 


3.3  Stability  Constraints 

There  are  two  time  steps  that  a  user  must  specify  when  running  the  ARW:  the  model  time  step 
(the  time  step  used  by  the  RK3  scheme,  see  Section  3.1.1)  and  the  acoustic  time  step  (used 
in  the  acoustic  sub-steps  of  the  time-split  integration  procedure,  see  Section  3.1.2).  Both  are 
limited  by  Courant  numbers.  In  the  following  sections  we  describe  how  to  choose  time  steps  for 
applications. 

3.3.1  RK3  Time  Step  Constraint 

The  RK3  time  step  is  limited  by  the  advective  Courant  number  uAt/Ax  and  the  user’s  choice 
of  advection  schemes —  users  can  choose  2”'’*  through  6*^  order  discretizations  for  the  advection 
terms.  The  time-step  limitations  for  ID  advection  in  the  RK3  scheme  using  these  advection 
schemes  is  given  in  Wicker  and  Skamarock  (2002),  and  is  reproduced  here. 


Time  Scheme 

Spatial  order 

3rd 

4th 

5th 

6th 

Leapfrog 

Unstable 

0.72 

Unstable 

0.62 

RK2 

0.88 

Unstable 

0.30 

Unstable 

RK3 

1.61 

1.26 

1.42 

1.08 

Table  3.1:  Maximum  stable  Courant  numbers  for  one-dimensional  linear  advection.  From  Wicker 
and  Skamarock  (2002). 

As  is  indicated  in  the  table,  the  maximum  stable  Courant  numbers  for  advection  in  the  RK3 
scheme  are  almost  a  factor  of  two  greater  than  those  for  the  leapfrog  time-integration  scheme. 
For  advection  in  three  spatial  dimensions,  the  maximum  stable  Courant  number  is  1/ a/S  times 
the  Courant  numbers  given  in  Table  3.1.  For  stability,  the  time  step  used  in  the  ARW  should 
produce  a  maximum  Courant  number  less  than  that  given  by  theory.  Thus,  for  3D  applications, 
the  time  step  should  satisfy  the  following  equation: 

(3.55) 

V  o  "^max 


21 


where  Cvtheory  is  the  Courant  number  taken  from  the  RK3  entry  in  Table  3.1  and  Umax  is  the 
maximum  velocity  expected  in  the  simulation.  For  example  in  real-data  applications,  where  jet 
stream  winds  may  reach  as  high  as  100  ms“^,  the  maximum  time  step  would  be  approximately 
80  s  on  a  Ax  =  10  km  grid  using  5*^  order  advection.  Given  additional  constraint  from  the  time 
splitting,  and  to  provide  a  safety  buffer,  we  usually  choose  a  time  step  that  is  approximately 
25%  less  than  that  given  by  (3.55).  This  time  step  is  typically  a  factor  of  two  greater  than  that 
used  in  leapfrog-based  models.  For  those  users  familiar  with  the  MM5  model,  the  rule  of  thumb 
for  choosing  a  time  step  is  that  the  time  step,  in  seconds,  should  be  approximately  3  times 
the  horizontal  grid  distance,  in  kilometers.  For  the  ARW,  the  time  step  (in  seconds)  should  be 
approximately  6  times  the  grid  distance  (in  kilometers). 

3.3.2  Acoustic  Time  Step  Constraint 

The  forward-backward  time  integration  scheme  used  in  the  ARW’s  2D  explicit  acoustic  step 
integration  allows  a  maximum  Courant  number  Cvmax  =  CgAr/Ax  <  l/\/2,  where  c*  is  the 
speed  of  sound.  We  typically  use  a  more  conservative  estimate  for  this  by  replacing  the  limiting 
value  l/-\/2  with  1/2.  Thus,  the  acoustic  time  step  used  in  the  model  is 

/\t 

At  <2-—.  (3.56) 

Cs 

For  example,  on  a  Ax  =  10  km  grid,  using  a  sound  speed  Cg  =  300  ms“^,  the  acoustic  time 
step  given  in  (3.56)  is  approximately  17  s.  In  the  ARW,  the  ratio  of  the  RK3  time  step  to  the 
acoustic  time  step  must  be  an  even  integer.  For  our  example  using  a  Ax  =  10  km  grid  in  a 
real-data  simulation,  we  would  specify  the  RK3  time  step  At  =  60s  (i.e.,  25%  less  than  the  80  s 
step  given  by  (3.55),  and  an  acoustic  time  step  At  =  15  s  (i.e.,  1/4  of  the  RK3  step,  rounding 
down  the  At  =  17  s  step  given  by  (3.56)).  Note  that  it  is  the  ratio  of  the  RK3  time  step  to  the 
acoustic  time  step  that  is  the  required  input  in  the  ARW. 
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Chapter  4 

Turbulent  Mixing  and  Model  Filters 


A  number  of  formulations  for  turbulent  mixing  and  filtering  are  available  in  the  ARW  solver. 
Some  of  these  filters  are  used  for  numerical  reasons.  For  example,  divergence  damping  is  used  to 
hlter  acoustic  modes  from  the  solution.  Other  filters  are  meant  to  represent  sub-grid  turbulence 
processes  that  cannot  be  resolved  on  the  chosen  grid.  These  filters  remove  energy  from  the 
solution  and  are  formulated  in  part  on  turbulence  theory  and  observations,  or  represent  energy 
sink  terms  in  some  approximation  to  the  Euler  equation.  In  this  section,  we  begin  by  outlining 
the  formulation  and  discretization  of  turbulent  mixing  processes  in  the  ARW  solver  commonly 
associated  with  sub-gridscale  turbulence  as  parameterized  in  cloud-scale  models —  the  second- 
order  horizontal  and  vertical  mixing.  In  large-scale  models  and  most  NWP  models,  vertical 
mixing  is  parameterized  within  the  planetary  boundary  layer  (PEL)  physics.  Vertical  mixing 
parameterized  within  the  PEL  physics  is  described  later  in  Chapter  8.  Here  we  note  that,  when 
a  PEL  parameterization  is  used,  all  other  vertical  mixing  is  disabled.  Following  the  outline  of 
turbulent  mixing  parameterizations  in  this  chapter,  other  numerical  hlters  available  in  the  ARW 
solver  are  described. 


4.1  Explicit  Spatial  Diffusion 

The  ARW  solver  has  two  formulations  for  spatial  dissipation —  diffusion  along  coordinate  sur¬ 
faces  and  diffusion  in  physical  (x,  y,  z)  space.  In  the  following  sections  we  present  the  diffusion 
operators  for  these  two  formulations,  followed  by  the  four  separate  formulations  that  can  be  used 
to  compute  the  eddy  viscosities.  We  conclude  with  a  description  of  the  prognostic  turbulent 
kinetic  energy  (TKE)  equation  used  in  one  set  of  these  formulations. 


4.1.1  Horizontal  and  Vertical  Diffusion  on  Coordinate  Surfaces 

For  any  model  variable,  horizontal  and  vertical  second  order  spatial  filtering  on  model  coordinate 
surfaces  is  considered  part  of  the  RHS  terms  in  the  continuous  equations  (2.35)  -  (2.41)  and  can 
be  expressed  as  follows  for  a  model  variable  a: 

dtifidd)  =  ...  +  fidirnd^imKhda^a)  +  mdy{mKhdya)]  +  {fid<y)~^  dy{K^a~^  d^a) .  (4.1) 
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For  the  horizontal  and  vertical  momentum  equations,  (4.1)  is  spatially  discretized  as 

dtU  =  ...  +  /I/m"  [6^{mKM  +  +  g\jl/a^)-^6^{K,{a^^)-^6r,u) 

dtV  =  ...  +  Jl/my  +  6y{mK,,6yv)]  + 

dtW  =  ...  +  iJ,dm[6a;im'^Kh''’6xw)  +  6y{m^  K  6yw)]  + 

The  spatial  discretization  for  a  scalar  q,  defined  at  the  mass  points,  is 
dtilJ^dq)  =  ■■■  +/idm[4(Fl"P“^iF/(5„g)  +  dyim^ P~^ K h  6yq)]  + 

In  the  current  ARW  formulation  for  mixing  on  coordinate  surfaces,  the  horizontal  eddy  viscosity 
Kh  is  allowed  to  vary  in  space,  whereas  the  vertical  eddy  viscosity  does  not  vary  in  space;  hence 
there  is  no  need  for  any  spatial  averaging  of  Additionally,  note  that  the  horizontal  eddy 
viscosity  Kh  is  multiplied  by  the  inverse  turbulent  Prandtl  number  P~^  for  horizontal  scalar 
mixing. 

4.1.2  Horizontal  and  Vertical  Diffusion  in  Physical  Space 

Coordinate  Metrics 

We  use  the  geometric  height  coordinate  in  this  physical  space  formulation.  The  coordinate 
metrics  are  computed  using  the  prognostic  geopotential  in  the  ARW  solver.  At  the  beginning 
of  each  Runge-Kutta  time  step,  the  coordinate  metrics  are  evaluated  as  part  of  the  overall 
algorithm.  The  definitions  of  the  metrics  are 

=  g~^Sx(j)  and  Zy  =  g~^6y(j). 

These  metric  terms  are  defined  on  w  levels,  and  {zx,Zy)  are  horizontally  coincident  with  {u,v) 
points.  Additionally,  the  vertical  diffusion  terms  are  evaluated  directly  in  terms  of  the  geometric 
height,  avoiding  the  need  for  metric  terms  in  the  vertical. 

Continuous  Equations 

The  continuous  equations  for  evaluating  diffusion  in  physical  space,  using  the  velocity  stress 
tensor,  are  as  follows  for  horizontal  and  vertical  momentum: 


dtU  =  .. 

.  -m[d^Tii  +  dyTu  -  d^{zxTii  +  %ri2)]  -  d^ris 

(4.2) 

dtV  =  .. 

.  -m[dxTi2  +  dyT22  -  dz{z^Ti2  +  ZyT22)]  “  dzT23 

(4.3) 

dtW  =  .. 

.  -m[d^Ti3  +  dyT23  -  dz{ZxTi3  +  ZyT23)]  - 

(4.4) 

The  stress  tensor  r  can  be  written  as  follows: 

Ri  =  —qi-dKhDii 

T22  =  —g-dKhD22 

T33  =  —g-dKjjD^s 
R2  =  —fJ'dKhDi2 
R3  =  —g^dKhDis 
R3  =  —fJ'dKhD23. 
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Symmetry  sets  the  remaining  tensor  values;  T21  =  T12,  =  Tis,  and  7^,2  =  T23.  The  stress  tensor 

r  is  calculated  from  the  deformation  tensor  D.  The  continuous  deformation  tensor  is  dehned  as 

Dll  =  ‘2rn?\dx{rn~^u)  —  Zxdz{rn~^u)\ 

D22  =  2m^[dy{m~^v)  -  Zydz{rrr^v)\ 

T>33  =  2  d^w 

Di2  =  \dy{m~^u)  —  Zydz{m~^u)  +  dx{rn~^v)  —  Zxdz{rn~^v)~\ 

Di3  =  [d^{m~^w)  -  z^d^{m~^w)]+d^{u) 

D23  =  rri^  [dy{m~^w)  -  Zydz{m~^w)]+d^{v). 

The  deformation  tensor  is  symmetric,  hence  D21  =  -D12,  -D31  =  D13,  and  D32  =  D23. 

The  diffusion  formulation  for  scalars  is 

dtifJ^dq)  =  •••  +  [m{d^  -  d^z^)  {ndmK{d^  -  z^d;.))  + 

m{dy  -  d^Zy)  {ndmK{dy  -  Zyd^))  +  d^/idKd^]  q.  (4.5) 


Spatial  Discretization 

Using  the  definition  of  the  stress  tensor,  the  spatial  discretization  of  the  ARW  physical-space 
diffusion  operators  for  the  horizontal  and  vertical  momentum  equations  (4.2)  -  (4.4)  are 


dtU  =  .. 

.  -  [dxTii  5yTi2  -  SziZxTll'"'’  Zy'"yTi2^'^)] 

—  13 

dtV  =  .. 

.  -  my  [6yr22  +  6xTi2  -  Sz^Zyf^yy  +  Zx^yiS^y)] 

—  ^^'r23 

dtW  =  .. 

.  -  m  [5xTi3  +  5yT23  - 

']  —  ^zT33 

The  discrete  forms  of  the  stress  tensor  and  deformation  tensor  are 


and 


'Til  —  —^idKhDii 

'r22  =  —fJ'dKhD22 


T33  —  —fJ'dKvD33 

ri2  = 

Ti3  =  -Wd^K~h’^Di3 
r23  =  -Jh^^''D23, 


Dll 

D22 

D33 

Di2 

Di3 

D23 


2m'^[6x{Tn  ^  u)  —  Zx^'^6z{m  ^  u)] 

2m^\5y{m~^^v)  — 'Zy^'^  5  v)\ 

2  6zW 


{m^yf 

6y{m  —  Zy^^dzim 

rn? 

m~^w)  —  ZxSz{m~^w)  ^ 

m? 

5y{m  ^w)  —  Zy5z{rn~^w)^^ 

_ 

-yv 


+6zU 

+5zV. 


Zx^'^5z{m  i^u) 


25 


The  spatial  discretization  for  the  scalar  diffusion  (4.5)  is 


dt{M)  =  ••• 

+  m[6y(jlj^ H2{q))  - 
+  Hd5z  {Kj^5^q) , 


where 


H,{q)=m^Kh{5,q-zJz{(f^)), 

H2{q)  =  myKh[5yq  -  Zy5^{q'^'^)) . 

4.1.3  Computation  of  the  Eddy  Viscosities 

There  are  four  options  for  determining  the  eddy  viscosities  and  in  the  ARW  solver. 


External  specification  of  Kh  and 

Constant  values  for  Kh  and  Ky  can  be  input  in  the  ARW  namelist. 


Kh  determined  from  the  horizontal  deformation 


The  horizontal  eddy  viscosity  Kh  can  be  determined  from  the  horizontal  deformation  using  a 
Smagorinsky  hrst-order  closure  approach.  In  this  formulation,  the  eddy  viscosity  is  dehned  and 
discretized  as 


Kh  =  C? 


0.25(T) 


11 


D 


22) 


+  Dl 


-xy 


12 


1 

2 


The  deformation  tensor  components  have  been  dehned  in  the  previous  section.  The  length  scale 
I  =  (AxAi/)^/^  and  Cg  is  a  constant  with  a  typical  value  Cg  =  0.25.  For  scalar  mixing,  the 
eddy  viscosity  is  divided  by  the  turbulent  Prandtl  number  Py  that  typically  has  a  value  of  1/3 
(Deardorff,  1972).  This  option  is  most  often  used  with  a  planetary  boundary  layer  scheme  that 
independently  handles  the  vertical  mixing. 


3D  Smagorinsky  Closure 

The  horizontal  and  vertical  eddy  viscosities  can  be  determined  using  a  3D  Smagorinsky  turbu¬ 
lence  closure.  This  closure  specihes  the  eddy  viscosities  as 


where 


Kh,v  =  cl  ll  y  max 


-lA^2^1/2 


0.,D^  -  {P-^N^) 


D=- 

2 


D?i + Di, + Di, + {d^2^y + {D^ry + 


-xr)\  2 


-y»?\2 


(4.6) 


and  N  is  the  Brunt- Vaisala  frequency;  the  computation  of  N,  including  moisture  effects,  is 
outlined  in  Section  4.1.4. 
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If  Ax  is  less  than  the  user-specihed  critical  length  scale  Icr,  then  the  length  scale  used  for 
calculating  both  and  in  (4.6)  is  Ih^y  =  {AxAyAzy^^  (and  Kh  =  =  K).  If  Ax  is 

greater  than  an  critical  length  scale  Icr,  then  the  length  scale  Ih  =  AxAy  in  the  calculation 
of  the  horizontal  eddy  viscosity  Kh  using  (4.6),  and  1^  =  Az  for  the  calculation  of  the  vertical 
eddy  viscosity  using  (4.6). 

Additionally,  the  eddy  viscosities  for  scalar  mixing  are  divided  by  the  turbulent  Prandtl 
number  Pr  =  1/3. 

Prognostic  TKE  Closure 

For  the  predicted  turbulent  kinetic  energy  option  (TKE;  see  section  4.1.4),  the  eddy  viscosities 
are  computed  using 

^h,v  Chlh,v  y/^, 

where  e  is  the  turbulent  kinetic  energy  (a  prognostic  variable  in  this  scheme),  Ck  is  a  constant 
(typically  0.15  <  Ck  <  0.25),  and  I  is  a  length  scale. 

If  the  grid  spacing  Ax  is  less  than  the  critical  length  scale  Icr,  then 

lh,v  =  min[(AxA2/A^)^/^,  0.76A/e/A]  for  >  0, 

lh,v  =  {AxAyAzY^^  for  <  0 

(see  section  4.1.4  for  the  calculation  of  A^).  Both  the  horizontal  and  vertical  eddy  viscosities 
are  multiplied  by  an  inverse  turbulent  Prandtl  number  P~^  =  1  +  2l/{AxAyAzY^^  for  scalar 
mixing.  In  this  case  {Icr  >  the  horizontal  and  vertical  eddy  viscosities  are  equivalent. 

If  the  grid  spacing  Ax  is  greater  than  the  critical  length  scale  Icr,  then  Ih  =  Y AxAy  for  the 
calculation  of  Kh-  For  calculating  K^, 

Ir  =  mm[Az,0.76^/e/N~\  for  A^  >  0, 

ly  =  Az  for  A^  <  0. 

The  eddy  viscosity  used  for  mixing  scalars  is  divided  by  a  turbulent  Prandtl  number  Pr-  The 
Prandtl  number  is  1/3  for  the  horizontal  eddy  viscosity  Kh,  and  Pp^  =  1  +  21/ Az  for  the  vertical 
eddy  viscosity  Ky. 

4.1.4  TKE  equation  for  the  1.5  Order  Turbulence  Closure 

The  prognostic  equation  governing  the  evolution  of  the  turbulent  kinetic  energy  e  is 

dt{+d^)  +  (V  ■  Ve)^  =  +d{  shear  production  +  buoyancy  +  dissipation  ).  (4.7) 

The  time  integration  and  the  transport  terms  in  (4.7)  are  integrated  in  the  same  manner  as  for 
other  scalars  (as  described  in  Chapter  3).  The  right-hand  side  source  and  sink  terms  for  e  are 
given  as  follows. 

Shear  Production 

The  shear  production  term  in  (4.7)  can  be  written  as 

shear  production  =  KhDl^  +  KhD^^  +  +  KhD\2  ^  +  KyDf^  ^  -|-  KyD^,^^ . 
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Buoyancy 

The  buoyancy  term  in  the  TKE  equation  (4.7)  is  written  as 

buoyancy  =  —KyN'^, 


where  the  Brunt- Vaisala  frequency  N  is  computed  using  either  the  formula  for  a  moist  saturated 
or  unsaturated  environment: 


N^  =  g 
N^  =  g 


'  dOe  dqy, 
dz  dz 


ide 

9  dz 


+  1.61^ 

oz 


dqw 

dz 


if  qv  >  qs  or  gc  >  0.01  g/Kg; 
if  q^  <  qs  or  q^  <  0.01  g/Kg. 


The  coefficient  A  is  dehned  as 


A 


1  l.eULg^ 

n-1  ^  fidT 

1  I  eL'^qv  ’ 


where  q^  represents  the  total  water  (vapor  -|-  all  liquid  species  -|-  all  ice  species),  L  is  the  latent 
heat  of  condensation  and  e  is  the  molecular  weight  of  water  over  the  molecular  weight  of  dry 
air.  6e  is  the  equivalent  potential  temperature  and  is  dehned  as 


Oe 


9  1  + 


eLq^s 

CpT 


5 


where  q^s  is  the  saturation  vapor  mixing  ratio. 


Dissipation 

If  Ax  is  less  than  the  critical  length  scale  Icr,  the  dissipation  term  in  (4.7)  is 


dissipation  = - j — , 

where 

As 

As  =  {AxAyAzY^^,  and 

I  =  min[{AxAyAzY^^,  0.76\/e/iV] . 


If  Ax  is  greater  than  the  critical  length  scale  I 

dissipation  =  - 


cr,  the  dissipation  term  in  (4.7)  is 

2V2  6^/2 

'Ts  T’ 


where 


kz 

1  kz/lo' 


Iq  =  min 


ab  Jq'  \fez  dz 


lo^  \/edz 

ab  =  0.2,  and  k  =  0.4  is  the  von  Karman  constant. 
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4.2  Filters  for  the  Time-split  RK3  scheme 

Three  filters  are  used  in  the  ARW  time-split  RK3  scheme  apart  from  those  in  the  model  physics: 
three-dimensional  divergence  damping  (an  acoustic  model  hlter);  an  external-mode  hlter  that 
damps  vertically-integrated  horizontal  divergence;  and  off-centering  of  the  vertically  implicit 
integration  of  the  vertical  momentum  equation  and  geopoential  equation.  Each  of  these  is 
described  in  the  following  sections. 

4.2.1  Divergence  Damping 

The  damping  of  the  full  mass  divergence  is  a  hlter  for  acoustic  modes  in  the  ARW  solver.  This 
3D  mass  divergence  damping  is  described  in  Skamarock  and  Klemp  (1992).  The  hltering  is 
accomplished  by  forward  weighting  the  pressure  update  in  the  acoustic  step  loop  described  in 
Section  3.1.3,  step  (6).  The  linearized  equation  of  state  (3.5)  is  used  to  diagnose  the  pressure  at 
the  new  time  r  after  the  U'\  V" ,  /i((,  and  0"  have  been  advanced.  Divergence  damping  consists 
of  using  a  modihed  pressure  in  the  computation  of  the  pressure  gradient  terms  in  the  horizontal 
momentum  equations  in  the  acoustic  steps  (Equations  (3.7)  and  (3.8)).  Denoting  the  updated 
value  as  the  modihed  pressure  p*’"  used  in  (3.7)  and  (3.8)  can  be  written  as 

where  is  the  damping  coefficient.  This  modihcation  is  equivalent  to  inserting  a  horizontal 
dihusion  term  into  the  equation  for  the  3D  mass  divergence,  hence  the  name  divergence  damping. 
A  divergence  damping  coefficient  =  0.1  is  typically  used  in  the  ARW  applications,  independent 
of  time  step  or  grid  size. 

4.2.2  External  Mode  Filter 

The  external  modes  in  the  solution  are  damped  by  hltering  the  vertically-integrated  horizontal 
divergence.  This  damping  is  accomplished  by  adding  a  term  to  the  horizontal  momentum 
equations.  The  additional  term  added  to  (3.7)  and  (3.8)  are 


6rU"  =  .. 

'7e^a;(^T— Arhrf) 

(4.9) 

6rV"  =  .. 

•  ^edyidT-Ari^d)- 

(4.10) 

The  quantity  dr-AriJ^  is  the  vertically- integrated  mass  divergence  (i.e.,  (3.20))  from  the  previous 
acoustic  step  (that  is  computed  using  the  time  r  values  of  U  and  R),  and  ye  is  the  external 
mode  damping  coefficient.  An  external  mode  damping  coefficient  ye  =  0.01  is  typically  used  in 
the  ARW  applications,  independent  of  time  step  or  grid  size. 

4.2.3  Semi-Implicit  Acoustic  Step  Off-centering 

Forward-in-time  weighting  of  the  vertically-implicit  acoustic-time-step  terms  damps  instabilities 
associated  vertically-propagating  sound  waves.  The  forward  weighting  also  damps  instabilities 
associated  with  sloping  mode  levels  and  horizontally  propagating  sound  waves  (see  Durran  and 
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Klemp,  1983;  Dudhia,  1995).  The  off-centering  is  accomplished  by  using  a  positive  (non-zero) 
coefficient  (3  (3.19)  in  the  acoustic  time-step  vertical  momentum  equation  (3.11)  and  geopotential 
equation  (3.12).  An  off-centering  coefficient  /3  =  0.1  is  typically  used  in  the  ARW  applications, 
independent  of  time  step  or  grid  size. 


4.3  Other  Damping 


4.3.1  Gravity  Wave  Absorbing  Layer 


A  gravity-wave  absorbing  layer  is  available  in  the  ARW  solver.  The  absorbing  layer  increases  the 
second-order  horizontal  and  vertical  eddy  viscosities  in  the  absorbing  layer  using  the  following 
formulation: 


Kdh  = 


At 


73  cos 


ztop  -  ^7r\ 

2)' 


and 


^dv 


Az^ 

^73  cos 


Ztop  - 

2/ 


Here  7^  is  a  user-specihed  damping  coefficient,  Ztop  is  the  height  of  the  model  top  for  a  particular 
grid  column,  Zd  is  the  depth  of  the  damping  layer  (from  the  model  top),  and  Kdh  and  Kdv  are 
the  horizontal  and  vertical  eddy  viscosities  for  the  wave  absorbing  layer.  If  other  spatial  hlters 
are  being  used,  then  the  eddy  viscosities  that  are  used  for  the  second-order  horizontal  and 
vertical  eddy  viscosities  are  the  maximum  of  {Kh-,  Kdh)  and  {Ky,  Kdv)-  The  effect  of  this  hlter 
on  gravity  waves  is  discussed  in  Klemp  and  Lilly  (1978),  where  guidance  on  the  choice  of  the 
damping  coefficeint  7^  can  also  be  found. 


4.3.2  Rayleigh  Damping  Layer 

A  Rayleigh  damping  layer  is  also  available  in  the  ARW  solver.  This  scheme  applies  a  tendency 
to  M,  V,  w,  and  6  to  gradually  relax  the  variable  back  to  a  predetermined  reference  state  value, 

t{z)  (u-u), 

t{z)  (v-v), 

t{z)w, 

t(z)  (0-e). 

Overbars  indicate  the  reference  state  helds,  which  are  functions  of  2;  only  and  are  typically 
dehned  as  the  initial  horizontally  homogeneous  helds  in  idealized  simulations.  The  reference 
state  vertical  velocity  is  assumed  to  be  zero.  The  variable  r  dehnes  the  vertical  structure  of  the 
damping  layer,  and  has  a  form  similar  to  the  Rayleigh  damper  in  Durran  and  Klemp  (1983), 


du 

m 

dv 

m 

dw 

de 


t{z) 


— 7r  sin^ 


0 


Ztop-z\ 

Zd  ) 
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for  2;  >  {ztop  -  Zd); 
otherwise. 


where  7^  is  a  user  specified  damping  coeficient,  Ztop  is  the  height  of  the  model  top  for  a  particular 
grid  column,  and  Zd  is  the  depth  of  the  damping  layer  (from  the  model  top). 

Because  the  model  surfaces  change  height  with  time  in  the  ARW  solver,  the  reference  state 
values  at  each  grid  point  need  to  be  recalculated  at  every  time  step.  Thus,  a  linear  interpolation 
scheme  is  used  to  calculate  updated  reference  state  values  based  on  the  height  of  the  model 
levels  at  each  time  step. 

The  effect  of  this  filter  on  gravity  waves  is  discussed  in  Klemp  and  Lilly  (1978),  where 
guidance  on  the  choice  of  the  damping  coefficeint  7^  can  also  be  found. 

4.3.3  Vertical- Velocity  Damping 

This  is  also  called  w-damping.  In  semi-operational  or  operational  NWP  applications,  the  model 
robustness  can  be  improved  by  detecting  locations  where  the  vertical  motion  approaches  the 
limiting  Courant  number  for  stability,  and  applying  a  Rayleigh  damping  term  in  the  vertical 
momentum  equation  to  stabilize  the  motion.  This  term  is  non-physical  and  should  only  be  used 
in  the  situation  where  many,  or  long,  model  runs  are  being  done,  and  there  is  no  option  for  a 
re-run  with  a  shorter  time-step  if  a  failure  occurs  due  to  an  excessively  strong  updraft.  This 
might  be  the  case,  for  example,  in  an  operational  setting  where  real-time  forecasts  have  to  be 
produced  on  time  to  be  useful.  However,  if  this  term  activates  frequently,  consideration  should 
be  given  to  reducing  the  model  time-step. 

The  term  is  calculated  from 

^  VLdt 


If  Cr  >  Cr/d,  then 

dtW  =  ...  -  fid  sign(IT)7^(Cr  -  Crp), 

where  7^  is  the  damping  coefficient  (typically  0.3  ms“^),  and  Crp  is  the  activation  Courant 
number  (typically  1.0).  The  ARW  outputs  the  location  of  this  damping  when  it  is  active. 
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Chapter  5 

Initial  Conditions 


The  ARW  may  be  run  with  initial  conditions  that  are  dehned  analytically  for  idealized  simu¬ 
lations,  or  it  may  be  run  using  interpolated  data  from  either  a  large-scale  analysis  or  forecast 
for  real-data  cases.  Both  2D  and  3D  tests  cases  for  idealized  simulations  are  provided.  Several 
sample  cases  for  real-data  simulations  are  provided,  which  rely  on  pre-processing  from  an  ex¬ 
ternal  package  that  converts  the  large-scale  GriB  data  into  a  format  suitable  for  ingest  by  the 
ARW’s  real-data  processor. 

The  programs  that  generate  the  specihc  initial  conditions  for  the  selected  idealized  or  real- 
data  case  function  similarly.  They  provide  the  ARW  with; 

•  input  data  that  is  on  the  correct  horizontal  and  vertical  staggering; 

•  hydrostatically  balanced  reference  state  and  perturbation  helds;  and 

•  metadata  specifying  such  information  as  the  date,  grid  physical  characteristics,  and  pro¬ 
jection  details. 

For  neither  the  idealized  nor  the  real-data  cases  are  the  initial  conditions  enhanced  with  obser¬ 
vations.  However,  output  from  the  ARW  system  initial  condition  programs  is  suitable  as  input 
to  the  WRF  variational  assimilation  package  (see  Chapter  9). 


5.1  Initialization  for  Idealized  Conditions 

The  ARW  comes  with  a  number  of  test  cases  using  idealized  environments,  including  moun¬ 
tain  waves  (em_hill2d_x),  squall  lines  (em_squall2d_x,  em_squall2d_y),  supercell  thunderstorms 
(em_quarter_ss),  gravity  currents  (em_grav2d_x),  and  baroclinic  waves  (em_b_wave).  A  brief  de¬ 
scription  of  these  test  cases  can  be  found  in  the  READMETestmases  hie  provided  in  the  ARW 
release.  The  test  cases  include  examples  of  atmospheric  hows  at  hne  scales  (e.g.,  the  gravity 
current  example  has  a  grid-spacing  of  100  meters  and  a  time  step  of  1  second)  and  examples  of 
how  at  large  scales  (e.g.,  the  baroclinic  wave  test  case  uses  a  grid-spacing  of  100  km  and  a  time 
step  of  600  s),  in  addition  to  the  traditional  mesoscale  and  cloudscale  model  simulations.  The 
test  suite  allows  an  ARW  user  to  easily  reproduce  these  known  solutions.  The  test  suite  is  also 
the  starting  point  for  constructing  idealized  how  simulations  by  modifying  initializations  that 
closely  resemble  a  desired  initialization. 

All  of  these  tests  use  as  input  a  ID  sounding  specihed  as  a  function  of  geometric  height 
(except  for  the  baroclinic  wave  case  that  uses  a  2D  sounding  specihed  in  [y^z]),  and,  with  the 
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exception  of  the  baroclinic  wave  test  case,  the  sounding  hies  are  in  text  format  that  can  be 
directly  edited  by  the  user.  The  ID  specihcation  of  the  sounding  in  these  test  hies  requires  the 
surface  values  of  pressure,  potential  temperature,  and  water  vapor  mixing  ratio,  followed  by  the 
potential  temperature,  vapor  mixing  ratio,  and  horizontal  wind  components  at  some  heights 
above  the  surface.  The  initialization  programs  for  each  case  assume  that  this  moist  sounding 
represents  an  atmosphere  in  hydrostatic  balance. 

Two  sets  of  thermodynamic  helds  are  needed  for  the  model —  the  reference  state  and  the 
perturbation  state  (see  Chapter  2  for  further  discussion  of  the  equations).  The  reference  state 
used  in  the  idealized  initializations  is  computed  using  the  input  sounding  from  which  the  mois¬ 
ture  is  discarded  (because  the  reference  state  is  dry).  The  perturbation  state  is  computed  using 
the  full  moist  input  sounding.  The  procedure  for  computing  the  hydrostatically-balanced  ARW 
reference  and  perturbation  state  variables  from  the  input  sounding  is  as  follows.  First,  density 
and  both  a  dry  and  full  hydrostatic  pressure  are  computed  from  the  input  sounding  at  the  input 
sounding  levels  z.  This  is  accomplished  by  integrating  the  hydrostatic  equation  vertically  up  the 
column  using  the  surface  pressure  and  potential  temperature  as  the  lower  boundary  condition. 
The  hydrostatic  equation  is 

SzP  =  +  iRd/Rv)qv''),  (5-1) 

where  is  a  two  point  average  between  input  sounding  levels,  and  6zP  is  the  difference  of 
the  pressure  between  input  sounding  levels  divided  by  the  height  difference.  Additionally,  the 
equation  of  state  is  needed  to  close  the  system: 


1 

Oid  —  —  — 
Pd 


Po 


Rd 

1  + 


(5.2) 


where  and  6  are  given  in  the  input  sounding.  (5.1)  and  (5.2)  are  a  coupled  set  of  nonlinear 
equations  for  p  and  p  in  the  vertical  integration,  and  they  are  solved  by  iteration.  The  dry 
pressure  on  input  sounding  levels  is  computed  by  integrating  the  hydrostatic  relation  down  from 
the  top,  excluding  the  vapor  component. 

Having  computed  the  full  pressure  (dry  plus  vapor)  and  dry  air  pressure  on  the  input  sound¬ 
ing  levels,  the  pressure  at  the  model  top  {pdht)  is  computed  by  linear  interpolation  in  height 
(or  possibly  extrapolation)  given  the  height  of  the  model  top  (an  input  variable).  The  column 
mass  pd  is  computed  by  interpolating  the  dry  air  pressure  to  the  surface  and  subtracting  it  from 
Pdht-  Given  the  column  mass,  the  dry-air  pressure  at  each  rj  level  is  known  from  the  coordinate 
dehnition  (2.1),  repeated  here 


P  =  {Pdh  -  Pdht)/ Pd  where  pd  =  Pdhs  -  Pdht, 


and  the  pressures  pdhs  and  pdht  refer  to  the  dry  atmosphere.  The  potential  temperature  from 
the  input  sounding  is  interpolated  to  each  of  the  model  pressure  levels,  and  the  equation  of  state 
(5.2)  is  used  to  compute  the  inverse  density  ad-  Finally,  the  ARW’s  hydrostatic  relation  (2.9), 
in  discrete  form 

^dPd 

is  used  to  compute  the  geopotential.  This  procedure  is  used  to  compute  the  reference  state 
(assuming  a  dry  atmosphere)  and  the  full  state  (using  the  full  moist  sounding).  The  perturbation 
variables  are  computed  as  the  difference  between  the  reference  and  full  state.  It  should  also  be 
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noted  that  in  the  nonhydrostatic  model  integration,  the  inverse  density  is  diagnosed  from 
the  geopotential  nsing  this  equation  of  state,  and  the  pressure  is  diagnosed  from  the  equation  of 
state  using  the  inverse  density  ad  and  the  prognostic  potential  temperature  9.  Thus,  the  ARW’s 
prognostic  variables  fid,  9,  and  0  are  in  exact  hydrostatic  balance  for  the  model  equations  (to 
machine  roundoff). 


5.2  Initialization  for  Real-Data  Conditions 

The  initial  conditions  for  the  real-data  cases  are  pre-processed  through  a  separate  package  called 
the  Standard  Initialization  (SI,  see  Fig.  5.1).  The  output  from  the  SI  is  passed  to  the  real- 
data  pre-processor  in  the  ARW —  program  real —  which  generates  initial  and  lateral  boundary 
conditions.  This  section  is  primarily  about  the  steps  taken  to  build  the  initial  and  the  lateral 
boundary  conditions  for  a  real-data  case.  Even  though  the  SI  is  outside  of  the  ARW  system,  a 
brief  description  is  appropriate  to  see  how  the  raw  meteorological  and  static  terrestrial  data  are 
brought  into  the  model  for  real-data  cases. 


5.2.1  Use  of  the  Standard  Initialization  by  the  ARW 

The  SI  is  a  set  of  programs  that  takes  terrestrial  and  meteorological  data  (typically  in  GriB 
format)  and  transforms  them  for  input  to  the  ARW  pre-processor  program  for  real-data  cases 
{real).  Figure  5.1  shows  the  flow  of  data  into  and  out  of  the  SI  system.  The  first  step  for  the 
SI  is  to  define  a  physical  grid  (including  the  projection  type,  location  on  the  globe,  number  of 
grid  points,  nest  locations,  and  grid  distances)  and  to  interpolate  static  fields  to  the  prescribed 
domain.  Independent  of  the  domain  configuration,  an  external  analysis  or  forecast  is  processed 
by  the  SFs  GriB  decoder,  which  diagnoses  required  fields  and  reformats  the  GriB  data  into  an 
internal  binary  format.  With  a  specified  domain,  the  SI  horizontally  interpolates  the  meteoro¬ 
logical  data  onto  the  projected  domain(s),  and  vertically  interpolates  the  data  to  the  ARW’s 
T]  coordinate.  The  output  data  from  the  SI  supplies  a  complete  3-dimensional  snapshot  of  the 
atmosphere  on  the  selected  model  grid’s  horizontal  and  vertical  staggering  at  the  selected  time 
slices,  which  is  sent  to  the  ARW  pre-processor  program  for  real-data  cases. 

The  input  to  the  ARW  real-data  processor  from  the  SI  contains  3-dimensionaI  fields  of  poten¬ 
tial  temperature  (K),  mixing  ratio  (kg/kg),  and  the  horizontal  components  of  momentum  (m/s, 
already  rotated  to  the  model  projection).  The  2-dimensional  static  terrestrial  fields  include: 
albedo,  Goriolis  parameters,  terrain  elevation,  vegetation/land-use  type,  land/water  mask,  map 
scale  factors,  map  rotation  angle,  soil  texture  category,  vegetation  greenness  fraction,  annual 
mean  temperature,  and  latitude/longitude.  The  2-dimensional  time-dependent  fields  from  the 
external  model,  after  processing  by  the  SI,  include:  fid  (Pa),  layers  of  soil  temperature  (K)  and 
soil  moisture  (kg/kg,  either  total  moisture,  or  binned  into  total  and  liquid  content),  snow  depth 
(m),  skin  temperature  (K),  and  fractional  sea  ice.  All  of  the  fields  in  the  final  output  from  the 
SI  are  on  the  correct  horizontal  and  vertical  staggering  for  the  ARW.  The  input  data  from  the 
SI  is  assumed  to  be  hydrostatically  balanced. 
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Figure  5.1:  Schematic  showing  the  data  flow  and  program  components  in  the  SI,  and  how  the 
SI  feeds  initial  data  to  the  ARW.  Letters  in  the  rectangular  boxes  indicate  program  names. 
GRID.GEN:  defines  the  model  domain  and  creates  static  flies  of  terrestrial  data.  GRIB_PREP: 
decodes  GriB  data.  HINTERP:  interpolates  meteorological  data  to  the  model  domain.  VIN- 
TERP:  vertically  interpolates  data  to  model  coordinate. 

5.2.2  Reference  and  Perturbation  State 

Identical  to  the  idealized  initializations,  there  is  a  partitioning  of  some  of  the  meteorological 
data  into  reference  and  perturbation  fields.  For  real-data  cases,  the  reference  state  is  defined  by 
terrain  elevation  and  the  following  three  constants: 

•  po  (10^  Pa)  reference  sea  level  pressure; 

•  To  (usually  270  to  300  K)  reference  sea  level  temperature;  and 

•  A  (50  K)  temperature  difference  between  the  pressure  levels  of  po  and  po/e. 

Using  these  parameters,  the  dry  reference  state  surface  pressure  is 

+  (5,3) 

From  (5.3),  the  3-dimensional  reference  pressure  is  computed  as  a  function  of  the  vertical  coor¬ 
dinate  r]  levels  and  the  model  top  pdht  (input  provided  by  SI  for  real-data  cases): 

Pd  =  'n  {Pdhs  -  Pdht)  +  Pdht-  (5.4) 

With  (5.4),  the  reference  temperature  is  a  straight  line  on  a  skew-T  plot,  defined  as 

T  =  To  +  A  In^. 

Po 
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From  the  reference  temperature  and  pressure,  the  reference  potential  temperature  is  then  dehned 
as 


?id 

Pd\  I  Po\^^ 


e,=  n  +  Ain^]  ^ 


Po/  \Pd. 

and  the  reciprocal  of  the  reference  density  using  (5.4)  and  (5.5)  is  given  by 


(5.5) 


_  1  RdOd  f Pd\ 

ttd  =  —  = 


,  .  (5-6) 

Pd  Po  \Po/ 

The  base  state  difference  of  the  dry  surface  pressure  from  (5.3)  and  the  model  top  is  given  as 

Pd  Pdhs  Pdht-  (^•'^) 

From  (5.6)  and  (5.7),  the  reference  state  geopotential  defined  from  the  hydrostatic  relation  is 

Pd" 


One  of  the  total  fields  provided  to  the  real-data  cases  by  the  SI  is  ^d-  The  perturbation  field 
given  the  reference  value  (5.7)  is 

Pd  =  Pd-  Pd-  (5.8) 

Starting  with  the  reference  state  helds  (5.4,  5.6,  and  5.7)  and  the  dry  surface  pressure  perturba¬ 
tion  (5.8),  the  perturbation  helds  for  pressure  and  inverse  density  are  diagnosed.  The  pressure 
perturbation  includes  moisture  and  is  diagnosed  from  the  hydrostatic  equation 


^vP  =  Pd  (^1  +  Pd^ 


which  is  integrated  down  from  at  the  model  top  (where  p'  =  0)  to  recover  p' .  The  total  dry 
inverse  density  is  given  as 


ad  =  —  9  1  -h  — g, 

Po  V  Rd 


P'd  +  Pd\ 
Po 


which  dehnes  the  perturbation  held  for  inverse  density 


^d  —  (^d  —  Old- 

The  perturbation  geopotential  is  diagnosed  from  the  hydrostatic  relation 

=  -  {PdOl'd  +  P'd^d) 

by  upward  integration  using  the  terrain  elevation  as  the  lower  boundary  condition.  In  future 
versions  of  the  real-data  pre-processor,  p'  will  be  re-diagnosed  consistent  with  the  method  used 
in  the  model  (2.21)  as  a  hnal  step.  No  modihcations  to  the  original  pd,  u,  v,  or  6  from  the 
SI  are  performed.  The  vertical  component  of  velocity  is  initialized  to  zero. 
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5.2.3  Generating  Lateral  Boundary  Data 

This  section  deals  with  generating  the  separate  lateral  boundary  condition  hie  used  exclusively 
for  the  real-data  cases.  For  information  on  which  lateral  boundary  options  are  available  for  both 
the  idealized  and  real-data  cases,  see  Chapter  (6). 

The  specihed  lateral  boundary  condition  for  the  coarse  grid  for  real-data  cases  is  supplied 
by  an  external  hie  that  is  generated  by  program  real.  This  hie  contains  records  for  the  helds  u, 
V,  9,  0',  and  that  are  used  by  the  ARW  to  constrain  the  lateral  boundaries  (other  helds 

are  in  the  boundary  hie,  but  are  not  used).  The  lateral  boundary  hie  holds  one  less  time  period 
than  was  processed  by  the  SI.  Each  of  these  variables  has  both  a  valid  value  at  the  initial  time 
of  the  lateral  boundary  time  and  a  tendency  term  to  get  to  the  next  boundary  time  period.  For 
example,  assuming  a  3-hourly  availability  of  data  from  the  SI,  the  hrst  time  period  of  the  lateral 
boundary  hie  for  u  would  contain  data  for  both  coupled  u  (map  scale  factor  and  fid  interpolated 
to  the  variable’s  staggering)  at  the  0  h  time 


and  a  tendency  value  dehned  as 


Uoh  = 


fn^ 


5 

O/i 


Ut  = 


u,h  -  U, 


Oh 


3h 


which  would  take  a  grid  point  from  the  initial  value  to  the  value  at  the  next  large-scale  time 
during  3  simulation  hours.  The  horizontal  momentum  helds  are  coupled  both  with  fid  and  the 
inverse  map  factor.  The  other  3-dimensional  helds  (6*,  qy,  and  0')  are  coupled  only  with  /i^.  The 
jj!^  lateral  boundary  held  is  not  coupled. 

Each  lateral  boundary  held  is  dehned  along  the  four  sides  of  the  rectangular  grid  (loosely 
referred  to  as  the  north,  south,  east,  and  west  sides).  The  boundary  values  and  tendencies  for 
vertical  velocity  and  the  non- vapor  moisture  species  are  included  in  the  external  lateral  boundary 
hie,  but  act  as  place-holders  for  the  nested  boundary  data  for  the  hne  grids.  The  width  of  the 
lateral  boundary  along  each  of  the  four  sides  is  user  selectable  at  run-time. 


5.2.4  Masking  of  Surface  Fields 

Some  of  the  meteorological  and  static  helds  are  “masked”.  A  masked  held  is  one  in  which  the 
values  are  typically  dehned  only  over  water  (e.g.,  sea  surface  temperature)  or  dehned  only  over 
land  (e.g.,  soil  temperature).  The  need  to  match  all  of  the  masked  helds  consistently  to  each 
other  requires  additional  steps  for  the  real-data  cases  due  to  the  masked  data’s  presumed  use  in 
various  physics  packages  in  the  soil,  at  the  surface,  and  in  the  boundary  layer.  If  the  land/water 
mask  for  a  location  is  hagged  as  a  water  point,  then  the  vegetation  and  soil  categories  must  also 
recognize  the  location  as  the  special  water  hag  for  each  of  their  respective  categorical  indices. 

The  values  for  the  soil  temperature  and  soil  moisture  come  from  the  SI  on  the  native  levels 
originally  dehned  for  those  variables  in  the  large-scale  model.  The  SI  does  no  vertical  inter¬ 
polation  for  the  soil  data.  While  it  is  typical  to  try  to  match  the  ARW  soil  scheme  with  the 
incoming  data,  that  is  not  a  requirement.  Pre-processor  real  will  vertically  interpolate  (linear  in 
depth  below  the  ground)  from  the  incoming  levels  to  the  requested  soil  layers  to  be  used  within 
the  model. 
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Chapter  6 

Lateral  Boundary  Conditions 


Several  lateral  boundary  condition  options  exist  for  the  ARW  that  are  suitable  for  idealized 
flows,  and  a  specihed  lateral  boundary  condition  for  real-data  simulations  is  available.  These 
choices  are  handled  via  a  run-time  option  in  the  Fortran  namelist  hie.  The  coarsest  grid  of 
any  single  simulation  is  eligible  for  any  of  the  lateral  boundary  selections.  For  example,  real- 
data  cases  could  use  combinations  of  periodic,  symmetric,  or  open  lateral  boundary  conditions 
instead  of  the  more  traditional  time- dependent  conditions  provided  by  an  external  boundary 
hie.  However,  use  of  the  specihed  time-dependent  lateral  boundary  conditions  for  one  of  the 
idealized  simulations  is  not  possible  because  an  external  boundary  hie  is  not  generated.  The 
ARW  supports  rectangular  horizontal  grid  rehnement  with  integer  ratios  of  the  parent  and  child 
grid  distances  and  time  steps. 

For  nesting,  all  hne  domains  use  the  nest  time-dependent  lateral  boundary  condition  where 
the  outer  row  and  column  of  the  hne  grid  is  specihed  from  the  parent  domain,  described  in 
Section  7.3. 


6.1  Periodic  Lateral  Boundary  Conditions 

Periodic  lateral  boundary  conditions  in  the  ARW  can  be  specihed  as  periodic  in  x  (west-east), 
y  (south- north) ,  or  doubly  periodic  in  {x,y).  The  periodic  boundary  conditions  constrain  the 
solutions  to  be  periodic;  that  is,  a  generic  model  state  variable  'ijj  will  follow  the  relation 

'ip{x  +  nL^,  y  +  mLy)  =  i){x,  y) 

for  all  integer  {n,m).  The  periodicity  lengths  (L^^Ly)  are 

[(dimension  of  the  domain  in  x)  -  l]Aa:  and  [(dimension  of  the  domain  in  y)  -  1]A?/. 

6.2  Open  Lateral  Boundary  Conditions 

Open  lateral  boundary  conditions,  also  referred  to  as  gravity-wave  radiating  boundary  condi¬ 
tions,  can  be  specihed  for  the  west,  east,  north,  or  south  boundary,  or  any  combination  thereof. 
The  gravity  wave  radiation  conditions  follow  the  approach  of  Klemp  and  Lilly  (1978)  and  Klemp 
and  Wilhelmson  (1978). 

There  are  a  number  of  changes  in  the  base  numerical  algorithm  described  in  Chapter  3  that 
accompany  the  imposition  of  these  conditions.  First,  for  the  normal  horizontal  velocities  along 
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a  boundary  on  which  the  condition  is  specihed,  the  momentum  equation  for  the  horizontal 
velocity,  (3.7)  or  (3.8),  is  replaced  by 


S^U"  =  -U*6^u, 

where  U*  =  min(7/  —  Cb,  0)  at  the  x  =  0  (western)  boundary,  U*  =  max(f/  +  Cb,  0)  at  the  x  =  L 
(eastern)  boundary,  and  likewise  for  the  south  and  north  boundaries  for  the  V  momentum.  The 
horizontal  difference  operator  6x  is  evaluated  in  a  one-sided  manner  using  the  difference  between 
the  boundary  value  and  the  value  one  grid-point  into  the  grid  from  the  boundary.  Cb  is  the  phase 
speed  of  the  gravity  waves  that  are  to  be  radiated;  it  is  specihed  as  a  model  constant  (for  more 
details  see  Klemp  and  Lilly,  1978;  Klemp  and  Wilhelmson,  1978). 

For  scalars  and  non-normal  momentum  variables,  the  boundary-perpendicular  hux  diver¬ 
gence  term  is  replaced  with 

6xiU^l^)  =  U*6x^l^  +  ^/J  6xU, 

where  U*  =  min(t/,  0)  at  the  x  =  0  -|-  Ax/2  (western)  scalar  boundary,  U*  =  max(17,  0)  at  the 
X  =  L  —  Ax/2  (eastern)  boundary,  and  likewise  for  the  south  and  north  boundaries  using  V.  As 
was  the  case  for  the  momentum  equations,  the  horizontal  difference  operator  6^  is  evaluated  in  a 
one-sided  manner  using  the  difference  between  the  boundary  value  and  the  value  one  grid-point 
into  the  grid  from  the  boundary. 

6.3  Symmetric  Lateral  Boundary  Conditions 

Symmetry  lateral  boundary  conditions  can  be  specihed  for  the  west,  east,  north,  or  south  bound¬ 
ary,  or  any  combination  thereof.  The  symmetry  boundaries  are  located  on  the  normal-velocity 
planes  at  the  lateral  edges  of  the  grids.  The  normal  velocities  are  zero  at  these  boundaries,  and 
on  either  side  of  the  boundary  the  normal  velocity  satishes  the  relation 

U±{xb  -  x)  =  -U±{xb  +  x), 

where  Xb  is  the  location  of  the  symmetry  boundary.  All  other  variables  satisfy  the  relation 

^jJ{Xb  -  x)  =  'lp{Xb  +  x). 


6.4  Specified  Lateral  Boundary  Conditions 

Primarily  for  real-data  cases,  the  specihed  boundary  condition  is  also  referred  to  as  a  relaxation, 
or  nudging,  boundary  condition.  There  are  two  uses  of  the  specihed  boundaries  in  the  ARW:  for 
the  outer-most  coarse  grid  or  for  the  time-dependent  boundaries  supplied  to  a  nested  grid.  The 
specihed  lateral  boundary  conditions  for  the  nest  are  automatically  selected  for  all  of  the  hne 
grids,  even  if  the  coarse  grid  is  using  combinations  of  the  symmetry,  periodic,  or  open  options. 
If  the  specihed  lateral  boundary  condition  is  selected  for  the  coarse  grid,  then  all  four  grid  sides 
(west,  east,  north,  and  south)  use  specihed  lateral  conditions. 

The  coarse  grid  specihed  lateral  boundary  is  comprised  of  both  a  specihed  and  a  relaxation 
zone  as  shown  in  Fig.  6.1).  For  the  coarse  grid,  the  specihed  zone  is  determined  entirely  by 
temporal  interpolation  from  an  external  forecast  or  analysis  (supplied  by  the  SI).  The  width 
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Real-Data  Lateral  Boundary  Condition:  Location  of  Specified  and  Relaxation  Zones 


East 


South 


Figure  6.1:  Specified  and  relaxation  zones  for  a  grid  with  a  single  specified  row  and  column, 
and  four  rows  and  columns  for  the  relaxation  zone.  These  are  typical  values  used  for  a  specified 
lateral  boundary  condition  for  a  real-data  case. 


of  the  specified  zone  is  run-time  configurable,  but  is  typically  set  to  1  (i.e.,  the  last  row  and 
column  along  the  outer  edge  of  the  coarse  grid).  The  second  region  of  the  lateral  boundary  for 
the  coarse  grid  is  the  relaxation  zone.  The  relaxation  zone  is  where  the  model  is  nudged  or 
relaxed  towards  the  large-scale  forecast  (e.g.,  rows  and  columns  2  through  5  in  Fig.  6.1).  The 
size  of  the  relaxation  zone  is  a  run-time  option. 

The  specihed  lateral  boundary  condition  for  the  coarse  grid  requires  an  external  file,  gener¬ 
ated  during  the  same  pre-processing  as  the  initial  condition  file.  Let  "0  be  any  prognostic  value 
having  a  lateral  boundary  entry,  after  Davies  and  Turner  (1977), 

=  Fi{'4)ls  -  V’)  -  F2A‘^{'iI)ls  -  V’),  (6.1) 

where  n  is  the  number  of  grid  points  in  from  the  outer  row  or  column  along  the  boundary 
{SpecZone  -|-  1  <  n  <  SpecZone  +  RelaxZone  —  1;  see  Fig.  6.1)  and  'ip^s  is  fhe  large-scale  value 
obtained  by  spatial  and  temporal  interpolation  from  the  analyses.  is  a  5-point  horizontal 
smoother  applied  along  r^-surfaces.  The  weighting  function  coefficients  Fi  and  F2  in  (6.1)  are 
given  by 

^  1  SpecZone  +  RelaxZone  —  n 

^  lOAt  RelaxZone  —  1  ’ 

^  1  SpecZone  +  RelaxZone  —  n 

^  50At  RelaxZone  —  1  ’ 

where  n  extends  only  through  the  relaxation  zone  {SpecZone+1  <  n  <  SpecZ eme+RelaxZ one— 
1).  Fi  and  F2  are  linear  ramping  functions  with  a  maximum  at  the  first  relaxation  row  or  column 
nearest  the  coarse  grid  boundary  (just  inside  of  the  specihed  zone). 
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On  the  coarse  grid,  the  specihed  boundary  condition  applies  to  the  horizontal  wind  com¬ 
ponents,  potential  temperature,  0',  and  water  vapor.  The  lateral  boundary  hie  contains 
enough  information  to  update  the  boundary  zone  values  through  the  entire  simulation  period. 
The  momentum  helds  are  coupled  with  and  the  inverse  map  factor  (both  at  the  specihc 
variable’s  horizontal  staggering  location),  and  the  other  3-dimensional  helds  are  coupled  with 
Hd-  The  jj!^  held  is  not  coupled  in  the  lateral  boundary  hie. 

Vertical  velocity  has  a  zero  gradient  boundary  condition  applied  in  the  specihed  zone  (usually 
the  outer-most  row  and  column).  Microphysical  variables,  except  vapor,  and  all  other  scalars 
have  how-dependent  boundary  conditions  applied  in  the  specihed  zone.  This  boundary  condition 
specihes  zero  on  inhow  and  zero-gradient  on  outhow.  Since  these  boundary  conditions  require 
only  information  from  the  interior  of  the  grid,  these  variables  are  not  in  the  specihed  boundary 
condition  hie. 
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Chapter  7 
Nesting 


The  ARW  supports  horizontal  nesting  that  allows  resolution  to  be  focused  over  a  region  of 
interest  by  introducing  an  additional  grid  (or  grids)  into  the  simulation.  In  the  current  im¬ 
plementation,  only  horizontal  rehnement  is  available:  there  is  no  vertical  nesting  option.  The 
nested  grids  are  rectangular  and  are  aligned  with  the  parent  (coarser)  grid  within  which  they  are 
nested.  Additionally,  the  nested  grids  allow  any  integer  spatial  {Axcoarse/ ‘^Xfine)  and  temporal 
refinements  of  the  parent  grid.  This  nesting  implementation  is  in  many  ways  similar  to  the 
implementations  in  other  mesoscale  and  cloudscale  models  (e.g.  MM5,  ARPS,  COAMPS).  The 
major  improvement  in  the  ARW’s  nesting  infrastruture  compared  with  techniques  used  in  other 
models  is  the  ability  to  compute  nested  simulations  efficiently  on  parallel  distributed-memory 
computer  systems,  which  includes  support  for  moving  nested  grids.  The  WRF  Software  Frame¬ 
work,  described  in  Michalakes  et  ah  (2004),  makes  these  advances  possible.  In  this  chapter  we 
describe  the  various  nesting  options  available  in  the  ARW  and  the  numerical  coupling  between 
the  grids. 

7.1  Overview 

l-Way  and  2- Way  Grid  Nesting 

Nested  grid  simulations  can  be  produced  using  either  1-way  nesting  or  2-way  nesting  as  outlined 
in  Fig.  7.1.  The  1-way  and  2- way  nesting  options  refer  to  how  a  coarse  grid  and  the  fine  grid 
interact.  In  both  the  1-way  and  2- way  simulation  modes,  the  fine  grid  boundary  conditions 
(i.e.,  the  lateral  boundaries)  are  interpolated  from  the  coarse  grid  forecast.  In  a  1-way  nest,  this 
is  the  only  information  exchange  between  the  grids  (from  coarse  grid  to  fine  grid).  Hence,  the 
name  1-way  nesting.  In  the  2- way  nest  integration,  the  fine  grid  solution  replaces  the  coarse  grid 
solution  for  coarse  grid  points  that  lie  inside  the  fine  grid.  This  information  exchange  between 
the  grids  is  now  in  both  directions  (coarse-to-fine  and  fine-to-coarse) .  Hence,  the  name  2-way 
nesting. 

The  1-way  nest  option  may  be  run  in  one  of  two  different  methods.  One  option  is  to 
produce  the  nested  simulation  as  two  separate  ARW  simulations  as  described  in  the  leftmost 
box  in  Fig.  7.1.  In  this  mode,  the  coarse  grid  is  integrated  first.  Output  from  the  coarse 
grid  integration  is  then  processed  to  provide  boundary  conditions  for  the  nested  run  (usually 
at  a  much  lower  temporal  frequency  than  the  coarse  grid  time  step),  and  this  is  followed  by 
the  complete  time  integration  of  fine  (nested)  grid.  Hence,  this  1-way  option  is  equivalent  to 
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Figure  7.1:  1-way  and  2-way  nesting  options  in  the  ARW. 


running  two  separate  simulations  with  a  processing  step  in  between.  Also  with  separate  grid 
simulations,  an  intermediate  re-analysis  (such  as  via  3D-Var,  see  Section  9)  can  be  used. 

The  second  1-way  option  (lockstep  with  no  feedback),  depicted  in  the  middle  box  in  Fig.  7.1, 
is  run  as  a  traditional  simulation  with  two  (or  more)  grids  integrating  concurrently,  except  with 
the  feedback  runtime  option  shut  off.  This  option  provides  lateral  boundary  conditions  to  the 
hne  grid  at  each  coarse  grid  time  step,  which  is  an  advantage  of  the  concurrent  1-way  method 
(no  feedback). 


Fine  Grid  Initialization  Options 

The  ARW  supports  several  strategies  to  refine  a  coarse-grid  simulation  with  the  introduction  of 
a  nested  grid.  When  using  1-way  and  2-way  nesting,  several  options  for  initializing  the  hne  grid 
are  provided. 

•  All  of  the  hne  grid  variables  can  be  interpolated  from  the  coarse  grid. 

•  All  of  the  hne  grid  variables  can  be  input  from  an  external  hie  which  has  high-resolution 
information  for  both  the  meteorological  and  the  terrestrial  helds. 

•  The  hne  grid  can  have  some  of  the  variables  initialized  with  a  high-resolution  external 
data  set,  while  other  variables  are  interpolated  from  the  coarse  grid. 

•  For  a  moving  nest,  an  external  orography  hie  may  be  used  to  update  the  hne  grid  terrain 
elevation.  This  option  is  not  generally  available  in  this  release. 

These  hne  grid  initialization  settings  are  user  specihed  at  run-time,  and  the  ARW  allows  nested 
grids  to  instantiate  and  cease  during  any  time  that  the  hne  grid’s  parent  is  still  integrating.  (The 
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Figure  7.2:  Various  nest  configurations  for  multiple  grids,  (a)  Telescoping  nests,  (b)  Nests  at 
the  same  level  with  respect  to  a  parent  grid,  (c)  Overlapping  grids:  not  allowed  (d)  Inner-most 
grid  has  more  than  one  parent  grid:  not  allowed 

system  is  currently  constrained  to  starting  nests  at  the  beginning  of  the  coarse  grid  simulations 
if  runs  require  input  of  nest-resolution  terrain  or  other  lower  boundary  data.  This  limitation 
will  be  addressed  in  the  near  future.) 

Possible  Grid  Configurations 

A  simulation  involves  one  outer  grid  and  may  contain  multiple  inner  nested  grids.  In  the  ARW, 
each  nested  region  is  entirely  contained  within  a  single  coarser  grid,  referred  to  as  the  parent 
grid.  The  hner,  nested  grids  are  referred  to  as  child  grids.  Using  this  terminology,  children 
are  also  parents  when  multiple  levels  of  nesting  are  used.  The  hne  grids  may  be  telescoped  to 
any  depth  (i.e.,  a  parent  grid  may  contain  one  or  more  child  grids,  each  of  which  in  turn  may 
successively  contain  one  or  more  child  grids;  Fig.  7.2a),  and  several  hne  grids  may  share  the 
same  parent  at  the  same  level  of  nesting  (Fig.  7.2b).  Any  valid  hne  grid  may  either  be  a  static 
domain  or  it  may  be  a  moving  nest  with  prescribed  incremental  shifts.  The  ARW  does  not 
permit  overlapping  grids,  where  a  coarse  grid  point  is  contained  within  more  than  a  single  child 
grid  (i.e.,  both  of  which  are  at  the  same  nest  level  with  respect  to  the  parent;  Fig.  7.2c).  In 
addition,  no  grid  can  have  more  than  a  single  parent  (Fig.  7.2d). 

For  2-way  nested  grid  simulations,  the  ratio  of  the  parent  horizontal  grid  distance  to  the 
child  horizontal  grid  distance  (the  spatial  rehnement  ratio)  must  be  an  integer.  This  is  also 
true  for  the  time  steps  (the  temporal  rehnement  ratio).  The  model  does  allow  the  time  step 
rehnement  ratio  to  diher  from  the  spatial  rehnement  ratio.  Also,  nested  grids  on  the  same  level 
(i.e.,  children  who  have  the  same  parent)  may  have  diherent  spatial  and  temporal  rehnement 
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Figure  7.3:  Arakawa-C  grid  staggering  for  a  portion  of  a  parent  domain  and  an  imbedded  nest 
domain  with  a  3:1  grid  size  ratio.  The  solid  lines  denote  coarse  grid  cell  boundaries,  and  the 
dashed  lines  are  the  boundaries  for  each  hne  grid  cell.  The  horizontal  components  of  velocity 
(“U”  and  “V”)  are  dehned  along  the  normal  cell  face,  and  the  thermodynamic  variables  (“0”) 
are  dehned  at  the  center  of  the  grid  cell  (each  square).  The  bold  typeface  variables  along  the 
interface  between  the  coarse  and  the  hne  grid  dehne  the  locations  where  the  specihed  lateral 
boundaries  for  the  nest  are  in  ehect. 


ratios. 


7.2  Nesting  and  Staggering 

The  ARW  uses  an  Arakawa-C  grid  staggering.  As  shown  in  Fig.  7.3,  the  u  and  v  com¬ 
ponents  of  horizontal  velocity  are  normal  to  the  respective  faces  of  the  grid  cell,  and  the 
mass/thermodynamic/scalar  variables  are  located  in  the  center  of  the  cell. 

The  variable  staggering  has  an  additional  column  of  u  in  the  x-direction  and  an  addi¬ 
tional  row  of  V  in  the  y-direction  because  the  normal  velocity  points  dehne  the  grid  bound¬ 
aries.  The  horizontal  momentum  components  rehect  an  average  across  each  cell-face,  while 
each  mass/thermodynamic/scalar  variable  is  the  representative  mean  value  throughout  the  cell. 
Feedback  is  handled  to  preserve  these  mean  values:  the  mass/thermodynamic/scalar  helds  are 
fed  back  with  an  average  from  within  the  entire  coarse  grid  point  (Fig.  7.3),  and  the  horizontal 
momentum  variables  are  averaged  along  their  respective  normal  coarse  grid  cell  faces. 

The  horizontal  interpolation  (to  instantiate  a  grid  and  to  provide  time-dependent  lateral 
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boundaries)  does  not  conserve  mass.  The  feedback  mechanism,  for  most  of  the  unmasked  fields, 
uses  cell  averages  (for  mass/thermodynamic/scalar  quantities)  and  cell-face  averages  for  the 
horizontal  momentum  fields. 

The  staggering  dehnes  the  way  that  the  hne  grid  is  situated  on  top  of  the  coarse  grid.  For  all 
odd  ratios  there  is  a  coincident  point  for  each  variable:  a  location  that  has  the  coarse  grid  and 
the  hne  grid  at  the  same  physical  point.  The  location  of  this  point  depends  on  the  variable.  In 
each  of  the  coarse-grid  cells  with  an  odd  ratio,  the  middle  hne-grid  cell  is  the  coincident  point 
with  the  coarse  grid  for  all  of  the  mass-staggered  helds  (Fig.  7.3).  For  the  horizontal  momentum 
variables  the  normal  velocity  has  coincident  points  along  the  grid  boundaries  for  odd  ratios. 

For  helds  that  are  averaged  back  to  the  coarse  grid  in  the  feedback,  the  mean  of  the  nine 
mass/thermodynamic/scalar  (for  example,  due  to  the  3:1  grid-distance  ratio  in  the  example 
shown  in  Fig.  7.3)  hne  grid  points  is  fed  back  to  the  coarse  grid.  These  helds  include  most 
3D  and  2D  arrays.  For  the  horizontal  momentum  helds  averaged  back  to  the  coarse  grid  in 
the  feedback,  the  mean  of  three  (for  example,  due  to  the  3:1  grid-distance  ratio  in  the  example 
shown  in  Fig.  7.3)  hne  grid  points  is  fed  back  to  the  coarse  grid  from  along  the  coincident  cell 
face.  The  helds  that  are  masked  due  to  the  land/sea  category  are  fed  back  directly  from  the 
coincident  points  for  odd  ratios.  Only  masked  helds  use  the  feedback  method  where  a  single 
point  from  the  hne  grid  is  assigned  to  the  coarse  grid. 

A  diherence  between  the  odd  and  even  grid-distance  ratios  is  in  the  feedback  from  the  hne 
grid  to  the  coarse  grid.  No  coincident  points  exist  for  the  single  point  feedback  mechanisms 
for  even  grid  distance  ratios  (such  as  used  for  the  land/sea  masked  2D  helds).  For  a  2:1  even 
grid  distance  ratio.  Figure  7.4  shows  that  each  coarse  grid  point  has  four  hne  grid  cells  that 
are  equally  close,  and  therefore  four  equally  eligible  grid  points  for  use  as  the  single  hne-grid 
point  that  feeds  back  to  the  coarse  grid.  The  single-point  feedback  is  arbitrarily  chosen  as  the 
south-west  corner  of  the  four  neighboring  points.  This  arbitrary  assignment  to  masked  helds 
implies  that  even  grid  distance  ratios  are  more  snited  for  idealized  simulations  where  masked 
helds  are  less  important. 

7.3  Nested  Lateral  Boundary  Conditions 

For  the  hne  grid  with  2-way  nesting  or  1-way  nesting  (using  a  concurrent  ARW  simulation, 
see  Fig.  7.1),  the  boundary  conditions  are  specihed  by  the  parent  grid  at  every  coarse-grid 
time  step.  The  nest  lateral  boundary  condition  behaves  similarly  to  the  specihed  boundary 
condition  for  real-data  cases  (see  Section  6.4),  but  the  relaxation  zone  is  not  active.  Prognostic 
variables  are  entirely  specihed  in  the  onter  row  and  column  of  the  hne  grid  throngh  spatial  and 
temporal  interpolation  from  the  coarse  grid  (the  coarse  grid  is  stepped  forward  in  time  prior  to 
advancement  of  any  child  grid  of  that  parent). 

7.4  Steps  to  Generate  a  Nest  Grid 

Only  the  concurrent  1-way  nest  option  or  the  2-way  nest  option  are  considered  in  this  section. 
The  1-way  nest  option  (using  two  consecutive  ARW  simulations,  see  Fig.  7.1)  is  functionally 
similar  to  two  separate,  single-grid  simulations  and  does  not  ht  the  following  description.  For  a 
multiple  grid  simulation  within  a  single  model  run,  there  are  some  additional  infrastructure  steps 
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that  are  required  (briefly  described  in  Fig.  7.5).  While  the  following  text  details  a  simulation 
with  a  single  coarse-grid  and  a  single  fine-grid,  this  implies  no  lack  of  generality  when  handling 
multiple  grid  levels  or  multiple  grids  on  the  same  level. 


Nest  Instantiation 

The  fine  grid  is  instantiated  as  a  child  of  a  parent  grid  at  the  requested  start  time.  This 
initialization  is  within  the  integration  step  for  the  parent  grid,  so  no  child  grid  can  begin  if 
the  parent  is  not  active.  To  All  in  the  correct  meteorological  fields,  an  initialization  routine  is 
called  to  horizontally  interpolate  the  coarse-grid  data  to  the  fine  grid  locations  using  a  monotone 
interpolation  scheme  (described  in  Smolarkiewicz  and  Grell,  1990)  for  most  fields  (i.e.,  the  same 
scheme  employed  for  generating  the  fine  grid  lateral  boundary  conditions)  and  a  simple  linear 
interpolation  or  averaging  scheme  for  masked  or  categorical  fields.  For  fields  that  are  masked 
with  the  land/sea  background  (such  as  land  only  fields  (e.g.,  snow),  or  water  only  fields  (e.g., 
sea  ice)),  the  interpolator  needs  to  know  what  field  defines  the  template  for  the  masking  (such 
as  the  land  use  category).  Part  of  the  automatic  code  generation  handles  calling  each  field  with 
its  associated  interpolator. 
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Integrate  Parent  Grid  One  Time  Step 

If  Nest  Grid  Start  Time 

(1)  Horizontally  Interpolate  Parent  to  Child  Grid 

(2)  Optionally  Input  High-Resolution  Child  Data 

(3)  Compute  Child  Reference  State 

(4)  Feedback  Child  Initial  Data  to  Parent  Grid 

(5)  Re-Compute  Parent  Reference  State 
End  If  Nest  Grid  Start  Time 

Solve  Time  Step  for  Parent  Grid  (see  Fig.  3.1) 

While  Existing  Nest  Grids  to  Integrate 

(1)  Lateral  Forcing  from  Parent  Grid  to  Child 

(2)  Integrate  Child  Grid  to  Current  Time  of  Parent  Grid 

(3)  Feedback  Child  Grid  Information  to  Parent  Grid 

End  While  Existing  Nest  Grids  to  Integrate 

End  Grid  Integrate 


Figure  7.5:  Nest  grid  integration  sequence. 


Fine  Grid  Input 

After  the  horizontal  interpolation  is  completed,  a  few  orographic-based  variables  are  saved  so 
that  they  may  be  used  to  blend  the  lateral  boundaries  along  the  coarse-grid/fine-grid  interface. 
The  terrain  elevation,  and  the  reference  geopotential  (0)  are  stored  for  later  use.  The  helds 
selected  as  input  from  the  hne  grid  input  hie  (for  the  concurrent  1-way  and  2-way  forecast 
methods  shown  in  Fig.  7.1)  are  ingested,  and  they  overwrite  the  arrays  that  were  horizontally 
interpolated  from  the  coarse  grid.  No  quality  control  for  data  consistency  is  performed  for  the 
hne  grid  input.  All  such  masked  checks  are  completed  by  the  ARW  real-data  pre-processor  real. 

Interface  Blended  Orography 

When  the  hne  grid  data  has  been  input,  the  previously-saved  orographic-based  helds  are  blended 
across  the  four  outer  rows  and  columns  of  the  hne  grid.  The  blending  is  a  simple  linear  weighting 
between  the  interpolated  coarse-grid  values  (the  saved  data)  and  the  hne  grid  values  from  the 
input  hie.  The  weighting  scheme  is  given  as: 

•  row/column  1:  100%  interpolated  coarse  grid,  0%  hne  grid, 

•  row/column  2:  75%  interpolated  coarse  grid,  25%  hne  grid, 

•  row/column  3:  50%  interpolated  coarse  grid,  50%  hne  grid, 

•  row/column  4:  25%  interpolated  coarse  grid,  75%  hne  grid,  and 

•  row/column  5:  0%  interpolated  coarse  grid,  100%  hne  grid, 

where  the  row  or  column  nearest  the  outer  edge  takes  precedence  in  ambiguous  corner  zones. 
The  blended  arrays  are  required  to  compute  the  reference  state  for  the  hne  grid.  The  hrst  row 
and  column  (100%  interpolated  from  the  coarse  grid)  ensures  that  the  reference  state  for  the 
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coarse  grid  and  fine  grid  is  consistent  along  the  fine  grid  boundary  interface.  The  blending  along 
the  inner  rows  and  columns  ramps  the  coarse  grid  reference  state  to  the  fine  grid  reference  state 
for  a  smooth  transition  between  the  grids. 

Feedback 

So  that  the  coarse  grid  and  the  fine  grid  are  consistent  at  coincident  points,  the  fine  grid  values 
are  fed  back  to  the  coarse  grid.  There  are  two  available  options  for  feedback:  either  the  mean  of 
all  fine  grid  cells  contained  within  each  coarse  grid  cell  is  fed  back  (or  cell  faces  in  the  case  of  the 
horizontal  momentum  fields),  or  a  single-point  feedback  is  selected  for  the  masked  or  categorical 
fields. 

Subsequent  to  the  feedback  step,  the  coarse  grid  may  be  optionally  smoothed  in  the  area  of 
the  fine  grid.  Two  smoothers  are  available:  a  5-point  1-2-1  smoother  and  a  smoother-desmoother 
with  a  similar  stencil  size.  Both  the  feedback  and  the  smoothers  are  run  one  row  and  column 
in  from  the  interface  row  and  column  of  the  coarse  grid  that  provides  the  lateral  boundary 
conditions  to  the  fine  grid. 

Reference  State 

The  initial  feedback  when  the  nest  is  instantiated  ensures  that  the  coarse  grid  is  consistent 
with  the  fine  grid,  particularly  with  regards  to  elevation  and  the  reference  state  fields  inside  the 
blended  region,  and  for  such  terrestrial  features  as  coasts,  lakes,  and  islands.  The  adjustment  of 
the  elevation  in  the  coarse  grid  forces  a  base  state  recalculation.  The  fine-grid  needs  an  initial 
base  state  calculation,  and  after  the  terrain  feedback,  the  coarse  grid  is  also  in  need  of  a  base 
state  recalculation. 

Note  that  with  the  horizontal  interpolation  of  the  coarse  grid  to  the  fine  grid  and  the  feedback 
of  the  fine  grid  to  the  coarse  grid,  the  coarse  grid  base  state  is  recomputed  even  without  a  separate 
fine-grid  initial  data  file,  since  the  coarse  grid  topography  is  adjusted. 

With  the  completed  base  state  computations,  the  routines  return  back  to  the  integration 
step  for  the  coarse  and  fine  grids.  The  fine  grid  data  is  now  properly  initialized  for  integration 
and  can  be  advanced  forward  a  time  step. 

Integration 

The  integration  by  grid  is  recursive.  At  the  end  of  each  grid’s  time  step,  a  check  is  made  to 
determine  if  a  child  grid  exists  for  that  parent  and  if  the  current  time  is  bracketed  by  the  child’s 
start/end  time.  This  is  shown  in  Fig.  7.5.  The  integration  process  for  the  nest  (step  2  under 
the  while  loop)  is  recursively  calling  the  top  step  in  the  overall  sequence  as  a  coarse  grid  itself. 
This  is  a  “depth  first”  traversal  of  the  tree  of  grids.  If  a  child  grid  does  exist,  that  child  grid  is 
integrated  up  through  the  current  time  of  the  parent  grid. 
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Chapter  8 
Physics 


This  chapter  outlines  the  physics  options  available  in  the  ARW.  The  WRF  physics  options 
fall  into  several  categories,  each  containing  several  options.  The  physics  categories  are  (1) 
microphysics,  (2)  cumulus  parameterization,  (3)  planetary  boundary  layer  (PEL),  (4)  land- 
surface  model,  and  (5)  radiation.  Diffusion,  which  may  also  be  considered  part  of  the  physics, 
is  described  in  Chapter  4. 

The  physics  section  is  insulated  from  the  rest  of  the  dynamics  solver  by  the  use  of  physics 
drivers.  These  are  between  solver-dependent  routines:  a  pre-physics  preparation  and  post¬ 
physics  modihcations  of  the  tendencies.  The  physics  preparation  involves  hlling  arrays  with 
physics-required  variables  that  include  the  temperature,  pressure,  heights,  layer  thicknesses, 
and  other  state  variables  in  MKS  units  at  half-level  grid  points  and  on  full  levels.  The  velocities 
are  also  de-staggered  so  that  the  physics  part  is  independent  of  the  dynamical  solver’s  velocity 
staggering.  Physics  packages  compute  tendencies  for  the  velocity  components  (un-staggered), 
potential  temperature,  and  moisture  fields.  The  solver- dependent  post-physics  step  will  re¬ 
stagger  these  tendencies  as  necessary,  couple  tendencies  with  coordinate  metrics,  and  convert  to 
variables  or  units  appropriate  to  the  dynamics  solver. 

In  the  hrst  Runge-Kutta  step,  prior  to  the  acoustic  steps  (see  Fig.  3.1,  step(l)),  tendencies 
are  computed  for  radiation,  surface,  PEL,  and  cumulus  physics.  These  tendencies  are  then  held 
hxed  through  the  Runge-Kutta  steps.  Microphysics  is  computed  after  the  last  Runge-Kutta 
step  (see  Fig.  3.1,  step  (9))  in  order  to  maintain  proper  saturation  conditions  at  the  end  of  the 
time- step. 

The  initialization  of  the  physics  is  called  prior  to  the  first  model  step.  This  initialization 
may  include  reading  in  data  files  for  physics  tables  or  calculating  look-up  tables  of  functions. 
Each  physics  module  includes  an  initialization  routine  for  this  purpose.  Often  physics  packages 
will  have  many  of  their  own  constants  that  should  also  be  included  in  their  own  module,  while 
common  physical  constants  are  passed  in  from  the  physics  drivers. 


8.1  Microphysics 

Microphysics  includes  explicitly  resolved  water  vapor,  cloud,  and  precipitation  processes.  The 
model  is  general  enough  to  accommodate  any  number  of  mass  mixing-ratio  variables,  and  other 
moments  such  as  number  concentrations.  Four-dimensional  arrays  with  three  spatial  indices  and 
one  species  index  are  use  to  carry  such  scalars.  Memory,  i.e.,  the  size  of  the  fourth  dimension 
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Table  8.1:  Microphysics  Options 


Scheme 

Number  of 
Variables 

Ice-Phase 

Processes 

Mixed-Phase 

Processes 

Kessler 

3 

N 

N 

Purdue  Lin 

6 

Y 

Y 

WSM3 

3 

Y 

N 

WSM5 

5 

Y 

N 

WSM6 

6 

Y 

Y 

Eta  GCP 

2 

Y 

Y 

Thompson 

7 

Y 

Y 

in  these  arrays,  is  allocated  depending  on  the  needs  of  the  scheme  chosen,  and  advection  of  the 
species  also  applies  to  all  those  reqnired  by  the  microphysics  option.  In  the  cnrrent  version  of 
the  ARW,  microphysics  is  carried  ont  at  the  end  of  the  time-step  as  an  adjnstment  process,  and 
so  does  not  provide  tendencies.  The  rationale  for  this  is  that  condensation  adjnstment  shonld 
be  at  the  end  of  the  time-step  to  gnarantee  that  the  final  satnration  balance  is  accnrate  for  the 
npdated  temperatnre  and  moistnre.  However,  it  is  also  important  to  have  the  latent  heating 
forcing  for  potential  temperature  during  the  dynamical  sub-steps,  and  this  is  done  by  saving  the 
microphysical  heating  as  an  approximation  for  the  next  time-step  as  described  in  Section  3.1.4. 

Currently,  the  sedimentation  process  is  accounted  for  in  the  microphysics,  and  a  smaller 
time  step  is  allowed  to  calculate  the  vertical  flux  of  precipitation  to  prevent  instability.  The 
saturation  adjustment  is  also  included  inside  the  microphysics.  In  the  future,  however,  it  might 
be  separated  into  an  individual  subroutine  to  enable  the  remaining  microphysics  to  be  called 
less  frequently  than  the  model’s  advection  step  for  efficiency. 

Table  8.1  shows  a  summary  of  the  options  indicating  the  number  of  moisture  variables,  and 
whether  ice-phase  and  mixed-phase  processes  are  included.  Mixed-phase  processes  are  those 
that  result  from  the  interaction  of  ice  and  water  particles,  such  as  riming  that  produces  graupel 
or  hail.  As  a  general  rule,  for  grid  sizes  less  than  10  km,  where  updrafts  may  be  resolved, 
mixed-phase  schemes  should  be  used,  particularly  in  convective  or  icing  situations.  For  coarser 
grids  the  added  expense  of  these  schemes  is  not  worth  it  because  riming  is  not  likely  to  be  well 
resolved. 


8.1.1  Kessler  scheme 

This  scheme  (Kessler,  1969),  which  was  taken  from  the  COMMAS  model  (Wicker  and  Wilhelm- 
son,  1995),  is  a  simple  warm  cloud  scheme  that  includes  water  vapor,  cloud  water,  and  rain.  The 
microphysical  processes  included  are:  the  production,  fall,  and  evaporation  of  rain;  the  accretion 
and  autoconversion  of  cloud  water;  and  the  production  of  cloud  water  from  condensation. 
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8.1.2  Purdue  Lin  scheme 


Six  classes  of  hydrometeors  are  included:  water  vapor,  cloud  water,  rain,  cloud  ice,  snow,  and 
graupel.  All  parameterization  production  terms  are  based  on  Lin  et  al.  (1983)  and  Rutledge 
and  Hobbs  (1984)  with  some  modihcations,  including  saturation  adjustment  following  Tao  et  ah 
(1989)  and  ice  sedimentation.  This  is  a  relatively  sophisticated  microphysics  scheme  in  WRF, 
and  it  is  more  suitable  for  use  in  research  studies.  The  scheme  is  taken  from  the  Purdue  cloud 
model,  and  the  details  can  be  found  in  Chen  and  Sun  (2002). 

8.1.3  WRF  Single-Moment  3-class  (WSM3)  scheme 

This  scheme  follows  Hong  et  al.  (2004)  including  ice  sedimentation  and  other  new  ice-phase 
parameterizations  revised  from  the  older  NCEP3  scheme  (Hong  et  al.,  1998)  that  was  in  WRF 
Version  1.  A  major  difference  from  other  schemes  is  that  a  diagnostic  relation  is  used  for 
ice  number  concentration  that  is  based  on  ice  mass  content  rather  than  temperature.  Three 
categories  of  hydrometers  are  included:  vapor,  cloud  water/ice,  and  rain/snow.  As  with  Dudhia 
(1989),  this  is  a  so-called  simple-ice  scheme  wherein  the  cloud  ice  and  cloud  water  are  counted 
as  the  same  category.  They  are  distinguished  by  temperature:  namely,  cloud  ice  can  only  exist 
when  the  temperature  is  less  than  or  equal  to  the  freezing  point;  otherwise,  cloud  water  can 
exist.  The  same  condition  is  applied  to  rain  and  snow.  Though  the  ice  phase  is  included,  it  is 
considered  efficient  enough  for  using  in  operational  models. 

8.1.4  WSM5  scheme 

This  scheme  is  similar  to  the  WSM3  simple  ice  scheme.  However,  vapor,  rain,  snow,  cloud  ice, 
and  cloud  water  are  held  in  hve  different  arrays.  Thus,  it  allows  supercooled  water  to  exist,  and 
a  gradual  melting  of  snow  as  it  falls  below  the  melting  layer.  Details  can  be  found  in  Hong  et 
al.  (2004).  It  replaces  WRF  Version  I’s  NCEP5  scheme  (Hong  et  ah,  1998). 

8.1.5  WSM6  scheme 

The  six-class  scheme  extends  the  WSM5  scheme  to  include  graupel  and  its  associated  processes. 
Many  of  these  processes  are  parameterized  similarly  to  Lin  et  al.  (1983),  but  there  are  differences 
for  the  accretion  calculation  and  in  some  other  parameters.  The  freezing/melting  processes  are 
computed  during  the  fall-term  sub-steps  to  increase  accuracy  in  the  vertical  heating  prohle  of 
these  processes.  The  order  of  the  processes  is  also  optimized  to  decrease  the  sensitivity  of  the 
scheme  to  the  time  step  of  the  model.  As  with  WSM3  and  WSM5,  saturation  adjustment  follows 
Dudhia  (1989)  and  Hong  et  al.  (1998)  in  separately  treating  ice  and  water  saturation  processes, 
rather  than  a  combined  saturation  such  as  the  Purdue  Lin  (above)  and  Goddard  (Tao  et  ah, 
1989)  schemes. 

8.1.6  Eta  Grid-scale  Cloud  and  Precipitation  (2001)  scheme 

This  is  also  known  as  EGCPOl  or  the  Eta  Ferrier  scheme.  The  scheme  predicts  changes  in 
water  vapor  and  condensate  in  the  forms  of  cloud  water,  rain,  cloud  ice,  and  precipitation  ice 
(snow/graupel/sleet).  The  individual  hydrometeor  helds  are  combined  into  total  condensate. 
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and  it  is  the  water  vapor  and  total  condensate  that  are  advected  in  the  model.  Local  storage  ar¬ 
rays  retain  hrst-guess  information  that  extract  contributions  of  cloud  water,  rain,  cloud  ice,  and 
precipitation  ice  of  variable  density  in  the  form  of  snow,  graupel,  or  sleet.  The  density  of  precip¬ 
itation  ice  is  estimated  from  a  local  array  that  stores  information  on  the  total  growth  of  ice  by 
vapor  deposition  and  accretion  of  liquid  water.  Sedimentation  is  treated  by  partitioning  the  time- 
averaged  flux  of  precipitation  into  a  grid  box  between  local  storage  in  the  box  and  fall  out  through 
the  bottom  of  the  box.  This  approach,  together  with  modifications  in  the  treatment  of  rapid 
microphysical  processes,  permits  large  time  steps  to  be  used  with  stable  results.  The  mean  size 
of  precipitation  ice  is  assumed  to  be  a  function  of  temperature  following  the  observational  results 
of  Ryan  (1996).  Mixed-phase  processes  are  now  considered  at  temperatures  warmer  than  -30°C 
(previously  -10°C),  whereas  ice  saturation  is  assumed  for  cloudy  conditions  at  colder  tempera¬ 
tures.  Further  description  of  the  scheme  can  be  found  in  Sec.  3.1  of  the  November  2001  Technical 
Procedures  Bulletin  (TPB)  at  http://www.emc.ncep.noaa.gov/mmb/mmbpll/etal2tpb/  and  on 
the  COMET  page  at  http://meted.ucar.edu/nwp/pcu2/etapcpl.htm. 

8.1.7  Thompson  et  al.  scheme 

The  Thompson  et  al.  (2004)  microphysical  parameterization  scheme  includes  improvements  to 
the  earlier  bulk  scheme  of  Reisner  et  al.  (1998)  and  has  been  extensively  tested  and  compared 
with  both  idealized  case  studies  and  documented  real  case  studies  of  mid-latitude  wintertime 
observations.  The  scheme  includes  six  classes  of  moisture  species  plus  number  concentration 
for  ice  as  prognostic  variables.  The  scheme  was  designed  to  improve  the  prediction  of  freezing 
drizzle  events  for  aircraft  safety  warnings.  Generally  microphysical  parameterizations  have  had 
problems  of  overpredicting  the  amount  of  snow  and  graupel  fields  and  under  predicting  the  ice 
in  outflow  regions  and  often  not  accurately  predicting  freezing  drizzle.  Key  improvements  are 
the  following: 

•  Primary  ice  nucleation  as  in  Cooper  (1986),  replaces  the  Fletcher  (1962)  curve. 

•  Auto-conversion  as  in  Walko  et  al.  (1995),  replaces  the  Kessler  (1969)  scheme. 

•  A  generalized  gamma  distribution  for  graupel  replaces  the  exponential  distribution. 

•  The  associated  intercept  parameter  depends  on  mixing  ratio  instead  of  remaining  constant. 

•  Riming  growth  of  snow  must  exceed  depositional  growth  of  snow  by  a  factor  of  3  before 
rimmed  snow  transfers  into  the  graupel  category. 

•  The  intercept  parameter  of  the  snow  size  distribution  depends  on  temperature. 

•  The  intercept  parameter  for  the  rain  size  distribution  depends  on  rain  mixing  ratio,  thereby 
simulating  the  fall  velocity  of  drizzle  drops  as  well  as  raindrops. 


8.2  Cumulus  parameterization 

These  schemes  are  responsible  for  the  sub-grid-scale  effects  of  convective  and/or  shallow  clouds. 
The  schemes  are  intended  to  represent  vertical  fluxes  due  to  unresolved  updrafts  and  downdrafts 
and  compensating  motion  outside  the  clouds.  They  operate  only  on  individual  columns  where 
the  scheme  is  triggered  and  provide  vertical  heating  and  moistening  profiles.  Some  schemes 
additionally  provide  cloud  and  precipitation  held  tendencies  in  the  column,  and  future  schemes 
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Table  8.2:  Cumulus  Parameterization  Options 


Scheme 

Gloud 

Detrainment 

Type  of  scheme 

Glosure 

Kain-Fritsch 

Y 

Mass  hux 

GAPE  removal 

Betts-Miller-Janjic 

N 

Adjustment 

Sounding  adjustment 

Grell-Devenyi 

Y 

Mass  hux 

Various 

may  provide  momentum  tendencies  due  to  convective  transport  of  momentum.  The  schemes  all 
provide  the  convective  component  of  surface  rainfall. 

Cumulus  parameterizations  are  theoretically  only  valid  for  coarser  grid  sizes,  (e.g.,  greater 
than  10  km),  where  they  are  necessary  to  properly  release  latent  heat  on  a  realistic  time  scale 
in  the  convective  columns.  While  the  assumptions  about  the  convective  eddies  being  entirely 
sub-grid-scale  break  down  for  finer  grid  sizes,  sometimes  these  schemes  have  been  found  to  be 
helpful  in  triggering  convection  in  5-10  km  grid  applications.  Generally,  they  should  not  be 
used  when  the  model  can  resolve  the  convective  eddies  itself  (e.g.,  <  5  km  grid). 

Table  8.2  summarizes  the  basic  characteristics  of  the  available  cumulus  parameterization 
options  in  the  ARW. 

8.2.1  Kain-Pritsch 

The  modified  version  of  the  Kain-Fritsch  scheme  (KF-Eta)  is  based  on  Kain  and  Fritsch  (1990) 
and  Kain  and  Fritsch  (1993),  but  has  been  modified  based  on  testing  within  the  Eta  model.  As 
with  the  original  KF  scheme,  it  utilizes  a  simple  cloud  model  with  moist  updrafts  and  downdrafts, 
including  the  effects  of  detrainment,  entrainment,  and  relatively  simple  microphysics.  It  differs 
from  the  original  KF  scheme  in  the  following  ways: 

•  A  minimum  entrainment  rate  is  imposed  to  suppress  widespread  convection  in  marginally 
unstable,  relatively  dry  environments. 

•  Shallow  (non  precipitating)  convection  is  allowed  for  any  updraft  that  does  not  reach 
minimum  cloud  depth  for  precipitating  clouds;  this  minimum  depth  varies  as  a  function 
of  cloud-base  temperature. 

•  The  entrainment  rate  is  allowed  to  vary  as  a  function  of  low-level  convergence. 

•  Downdraft  changes: 

—  Source  layer  is  the  entire  150  -  200  mb  deep  layer  just  above  cloud  base. 

—  Mass  flux  is  specihed  as  a  fraction  of  updraft  mass  flux  at  cloud  base.  Fraction  is 
a  function  of  source  layer  RH  rather  than  wind  shear  or  other  parameters,  i.e.,  old 
precipitation  efficiency  relationship  not  used. 

—  Detrainment  is  specihed  to  occur  in  updraft  source  layer  and  below. 
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8.2.2  Betts-Miller-Janjic 

The  Betts-Miller-Janjic  (BMJ)  scheme  (Janjic,  1994,  2000)  was  derived  from  the  Betts-Miller 
(BM)  convective  adjustment  scheme  (Betts,  1986;  Betts  and  Miller,  1986).  However,  the  BMJ 
scheme  differs  from  the  Betts-Miller  scheme  in  several  important  aspects.  The  deep  convec¬ 
tion  prohles  and  the  relaxation  time  are  variable  and  depend  on  the  cloud  efficiency,  a  non- 
dimensional  parameter  that  characterizes  the  convective  regime  (Janjic,  1994).  The  cloud  effi¬ 
ciency  depends  on  the  entropy  change,  precipitation,  and  mean  temperature  of  the  cloud.  The 
shallow  convection  moisture  prohle  is  derived  from  the  requirement  that  the  entropy  change 
be  small  and  nonnegative  (Janjic,  1994).  The  BMJ  scheme  has  been  optimized  over  years  of 
operational  application  at  NCEP,  so  that,  in  addition  to  the  described  conceptual  differences, 
many  details  and/or  parameter  values  differ  from  those  recommended  in  Betts  (1986)  and  Betts 
and  Miller  (1986).  Recently,  attempts  have  been  made  to  rehne  the  scheme  for  higher  horizontal 
resolutions,  primarily  through  modihcations  of  the  triggering  mechanism.  In  particular: 

•  A  floor  value  for  the  entropy  change  in  the  cloud  is  set  up  below  which  the  deep  convection 
is  not  triggered; 

•  In  searching  for  the  cloud  top,  the  ascending  particle  mixes  with  the  environment;  and 

•  The  work  of  the  buoyancy  force  on  the  ascending  particle  is  required  to  exceed  a  prescribed 
positive  threshold. 

8.2.3  Grell-Devenyi  ensemble 

Grell  and  Devenyi  (2002)  introduced  an  ensemble  cumulus  scheme  in  which  effectively  multiple 
cumulus  schemes  and  variants  are  run  within  each  grid  box  and  then  the  results  are  averaged 
to  give  the  feedback  to  the  model.  In  principle,  the  averaging  can  be  weighted  to  optimize  the 
scheme,  but  the  default  is  an  equal  weight.  The  schemes  are  all  mass-flux  type  schemes,  but  with 
differing  updraft  and  downdraft  entrainment  and  detrainment  parameters,  and  precipitation 
efficiencies.  These  differences  in  static  control  are  combined  with  differences  in  dynamic  control, 
which  is  the  method  of  determining  cloud  mass  flux.  The  dynamic  control  closures  are  based  on 
convective  available  potential  energy  (CAPE  or  cloud  work  function),  low-level  vertical  velocity, 
or  moisture  convergence.  Those  based  on  CAPE  either  balance  the  rate  of  change  of  CAPE  or 
relax  the  CAPE  to  a  climatological  value,  or  remove  the  CAPE  in  a  convective  time  scale.  The 
moisture  convergence  closure  balances  the  cloud  rainfall  to  the  integrated  vertical  advection 
of  moisture.  Another  control  is  the  trigger,  where  the  maximum  cap  strength  that  permits 
convection  can  be  varied.  These  controls  typically  provide  ensembles  of  144  members. 

8.3  Surface  Layer 

The  surface  layer  schemes  calculate  friction  velocities  and  exchange  coefficients  that  enable 
the  calculation  of  surface  heat  and  moisture  fluxes  by  the  land-surface  models  and  surface 
stress  in  the  planetary  boundary  layer  scheme.  Over  water  surfaces,  the  surface  fluxes  and 
surface  diagnostic  helds  are  computed  in  the  surface  layer  scheme  itself.  The  schemes  provide  no 
tendencies,  only  the  stability-dependent  information  about  the  surface  layer  for  the  land-surface 
and  PBL  schemes.  Currently,  each  surface  layer  option  is  tied  to  particular  boundary-layer 
options,  but  in  the  future  more  interchangeability  and  options  may  become  available. 
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8.3.1  Similarity  theory  (MM5) 

This  scheme  uses  stability  functions  from  Paulson  (1970),  Dyer  and  Hicks  (1970),  and  Webb 
(1970)  to  compute  surface  exchange  coefficients  for  heat,  moisture,  and  momentum.  A  convective 
velocity  following  Beljaars  (1994)  is  used  to  enhance  surface  fluxes  of  heat  and  moisture.  No 
thermal  roughness  length  parameterization  is  included  in  the  current  version  of  this  scheme. 
A  Charnock  relation  relates  roughness  length  to  friction  velocity  over  water.  There  are  four 
stability  regimes  following  Zhang  and  Anthes  (1982).  This  surface  layer  scheme  must  be  run  in 
conjunction  with  the  MRF  or  YSU  PEL  schemes. 

8.3.2  Similarity  theory  (Eta) 

The  Eta  surface  layer  scheme  (Janjic,  1996,  2002)  is  based  on  similarity  theory  (Monin  and 
Obukhov,  1954).  The  scheme  includes  parameterizations  of  a  viscous  sub-layer.  Over  water 
surfaces,  the  viscous  sub-layer  is  parameterized  explicitly  following  Janjic  (1994).  Over  land, 
the  effects  of  the  viscous  sub-layer  are  taken  into  account  through  variable  roughness  height  for 
temperature  and  humidity  as  proposed  by  Zilitinkevich  (1995).  The  Beljaars  (1994)  correction 
is  applied  in  order  to  avoid  singularities  in  the  case  of  an  unstable  surface  layer  and  vanishing 
wind  speed.  The  surface  fluxes  are  computed  by  an  iterative  method.  This  surface  layer  scheme 
must  be  run  in  conjunction  with  the  Eta  (Mellor-Yamada-Janjic)  PBL  scheme,  and  is  therefore 
sometimes  referred  to  as  the  MYJ  surface  scheme. 


8.4  Land-Surface  Model 

The  land-surface  models  (LSMs)  use  atmospheric  information  from  the  surface  layer  scheme, 
radiative  forcing  from  the  radiation  scheme,  and  precipitation  forcing  from  the  microphysics 
and  convective  schemes,  together  with  internal  information  on  the  land’s  state  variables  and 
land-surface  properties,  to  provide  heat  and  moisture  fluxes  over  land  points  and  sea-ice  points. 
These  fluxes  provide  a  lower  boundary  condition  for  the  vertical  transport  done  in  the  PBL 
schemes  (or  the  vertical  diffusion  scheme  in  the  case  where  a  PBL  scheme  is  not  run,  such  as  in 
large-eddy  mode).  [Note  that  large-eddy  mode  with  interactive  surface  fluxes  is  not  yet  available 
in  the  ARW,  but  is  planned  for  the  near  future.]  The  land-surface  models  have  various  degrees  of 
sophistication  in  dealing  with  thermal  and  moisture  fluxes  in  multiple  layers  of  the  soil  and  also 
may  handle  vegetation,  root,  and  canopy  effects  and  surface  snow-cover  prediction.  The  land- 
surface  model  provides  no  tendencies,  but  does  update  the  land’s  state  variables  which  include 
the  ground  (skin)  temperature,  soil  temperature  profile,  soil  moisture  profile,  snow  cover,  and 
possibly  canopy  properties.  There  is  no  horizontal  interaction  between  neighboring  points  in  the 
LSM,  so  it  can  be  regarded  as  a  one-dimensional  column  model  for  each  WRF  land  grid-point, 
and  many  LSMs  can  be  run  in  a  stand-alone  mode.  Table  8.3  summarizes  the  basic  features  of 
the  land-surface  treatments  in  ARW. 

8.4.1  5-layer  thermal  diffusion 

This  simple  LSM  is  based  on  the  MM5  5-layer  soil  temperature  model.  Layers  are  1,  2,  4,  8,  and 
16  cm  thick.  Below  these  layers,  the  temperature  is  hxed  at  a  deep-layer  average.  The  energy 
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Table  8.3:  Land  Surface  Options 


Scheme 

Vegetation 

Processes 

Soil 

Variables  (Layers) 

Snow 

Scheme 

5- layer 

N 

Temperature  (5) 

none 

Noah 

Y 

Temperature,  Water -|-Ice,  Water  (4) 

1-layer,  fractional 

RUC 

Y 

Temperature,  Ice,  Water  -|-  Ice  (6) 

multi-layer 

budget  includes  radiation,  sensible,  and  latent  heat  flux.  It  also  allows  for  a  snow-cover  flag,  but 
the  snow  cover  is  hxed  in  time.  Soil  moisture  is  also  hxed  with  a  landuse-  and  season-dependent 
constant  value,  and  there  are  no  explicit  vegetation  effects. 

8.4.2  Noah  LSM 

The  Noah  LSM  is  the  successor  to  the  OSU  LSM  described  by  Chen  and  Dudhia  (2001).  The 
scheme  was  developed  jointly  by  NCAR  and  NCEP,  and  is  a  unihed  code  for  research  and 
operational  purposes,  being  almost  identical  to  the  code  used  in  the  NCEP  North  American 
Mesoscale  Model  (NAM).  This  has  the  beneht  of  being  consistent  with  the  time-dependent  soil 
helds  provided  in  the  analysis  datasets.  This  is  a  4-layer  soil  temperature  and  moisture  model 
with  canopy  moisture  and  snow  cover  prediction.  It  includes  root  zone,  evapotranspiration,  soil 
drainage,  and  runoff,  taking  into  account  vegetation  categories,  monthly  vegetation  fraction, 
and  soil  texture.  The  scheme  provides  sensible  and  latent  heat  fluxes  to  the  boundary-layer 
scheme.  The  Noah  LSM  additionally  predicts  soil  ice,  and  fractional  snow  cover  effects,  has  an 
improved  urban  treatment,  and  considers  surface  emissivity  properties,  which  are  all  new  since 
the  OSU  scheme. 

8.4.3  Rapid  Update  Cycle  (RUC)  Model  LSM 

This  is  a  LSM  with  6  sub-soil  layers  and  up  to  two  snow  layers  that  is  used  operationally  in  the 
RUC  model  (Smirnova  et  al.,  1997,  2000).  The  model  considers  frozen  soil  processes,  patchy 
snow,  with  snow  temperature  and  density  variation,  vegetation  effects,  and  canopy  water. 


8.5  Planetary  Boundary  Layer 

The  planetary  boundary  layer  (PEL)  is  responsible  for  vertical  sub-grid-scale  fluxes  due  to  eddy 
transports  in  the  whole  atmospheric  column,  not  just  the  boundary  layer.  Thus,  when  a  PEL 
scheme  is  activated,  explicit  vertical  diffusion  is  de-activated  with  the  assumption  that  the  PEL 
scheme  will  handle  this  process.  The  most  appropriate  horizontal  diffusion  choices  (Section 
4.1.3)  are  those  based  on  horizontal  deformation  or  constant  K^,  values  where  horizontal  and 
vertical  mixing  are  treated  independently.  The  surface  fluxes  are  provided  by  the  surface  layer 
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Table  8.4:  Planetary  Boundary  Layer  Options 


Scheme 

Unstable  PBL 

Mixing 

Entrainment 

treatment 

PBL  Top 

MRF 

K  prohle  -|-  countergradient  term 

part  of  PBL  mixing 

from  critical  bulk  Ri 

YSU 

K  prohle  -|-  countergradient  term 

explicit  term 

from  buoyancy  prohle 

MYJ 

K  from  prognostic  TKE 

part  of  PBL  mixing 

from  TKE 

and  land-surface  schemes.  The  PBL  schemes  determine  the  flux  prohles  within  the  well-mixed 
boundary  layer  and  the  stable  layer,  and  thus  provide  atmospheric  tendencies  of  temperature, 
moisture  (including  clouds),  and  horizontal  momentum  in  the  entire  atmospheric  column.  Most 
PBL  schemes  consider  dry  mixing,  but  can  also  include  saturation  effects  in  the  vertical  stability 
that  determines  the  mixing.  The  schemes  are  one-dimensional,  and  assume  that  there  is  a  clear 
scale  separation  between  sub-grid  eddies  and  resolved  eddies.  This  assumption  will  become  less 
clear  at  grid  sizes  below  a  few  hundred  meters,  where  boundary  layer  eddies  may  start  to  be 
resolved,  and  in  these  situations  the  scheme  should  be  replaced  by  a  fully  three-dimensional 
local  sub-grid  turbulence  scheme  such  as  the  TKE  diffusion  scheme  (Section  4.1.4.)  Table  8.4 
summarizes  the  basic  features  of  the  PBL  schemes  in  ARW. 

8.5.1  Medium  Range  Forecast  Model  (MRF)  PBL 

The  scheme  is  described  by  Hong  and  Pan  (1996).  This  PBL  scheme  employs  a  so-called 
counter-gradient  flux  for  heat  and  moisture  in  unstable  conditions.  It  uses  enhanced  vertical 
flux  coefficients  in  the  PBL,  and  the  PBL  height  is  determined  from  a  critical  bulk  Richardson 
number.  It  handles  vertical  diffusion  with  an  implicit  local  scheme,  and  it  is  based  on  local  Ri 
in  the  free  atmosphere. 

8.5.2  Yonsei  University  (YSU)  PBL 

The  Yonsei  University  PBL  is  the  next  generation  of  the  MRF  PBL,  also  using  the  counter¬ 
gradient  terms  to  represent  fluxes  due  to  non-local  gradients.  This  adds  to  the  MRF  PBL 
an  explicit  treatment  of  the  entrainment  layer  at  the  PBL  top.  The  entrainment  is  made 
proportional  to  the  surface  buoyancy  flux  in  line  with  results  from  studies  with  large-eddy 
models.  The  PBL  top  is  dehned  using  a  critical  bulk  Richardson  number  of  zero  (compared  to 
0.5  in  the  MRF  PBL),  so  is  effectively  only  dependent  on  the  buoyancy  prohle  which,  in  general, 
lowers  the  calculated  PBL  top  compared  to  MRF. 

8.5.3  Mellor-Yamada-Janjic  (MYJ)  PBL 

This  parameterization  of  turbulence  in  the  PBL  and  in  the  free  atmosphere  (Janjic,  1990,  1996, 
2002)  represents  a  nonsingular  implementation  of  the  Mellor-Yamada  Level  2.5  turbulence  clo- 
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Table  8.5:  Radiation  Options 


Scheme 

Longwave/ 

Shortwave 

Spectral 

Bands 

GO2,  O3,  clouds 

RRTM 

LW 

16 

CC)25  Os,  clouds 

GFDL  LW 

LW 

14 

CO2,  O3,  clouds 

GFDL  SW 

SW 

12 

CO2,  O3,  clouds 

MM5  SW 

SW 

1 

clouds 

Goddard 

SW 

11 

CO2,  O3,  clouds 

snre  model  (Mellor  and  Yamada,  1982)  through  the  full  range  of  atmospheric  turbulent  regimes. 
In  this  implementation,  an  upper  limit  is  imposed  on  the  master  length  scale.  This  upper  limit 
depends  on  the  TKE  as  well  as  the  buoyancy  and  shear  of  the  driving  flow.  In  the  unstable 
range,  the  functional  form  of  the  upper  limit  is  derived  from  the  requirement  that  the  TKE  pro¬ 
duction  be  nonsingular  in  the  case  of  growing  turbulence.  In  the  stable  range,  the  upper  limit 
is  derived  from  the  requirement  that  the  ratio  of  the  variance  of  the  vertical  velocity  deviation 
and  TKE  cannot  be  smaller  than  that  corresponding  to  the  regime  of  vanishing  turbulence.  The 
TKE  production/dissipation  differential  equation  is  solved  iteratively.  The  empirical  constants 
have  been  revised  as  well  (Janjic,  1996,  2002). 


8.6  Atmospheric  Radiation 


The  radiation  schemes  provide  atmospheric  heating  due  to  radiative  flux  divergence  and  surface 
downward  longwave  and  shortwave  radiation  for  the  ground  heat  budget.  Longwave  radiation 
includes  infrared  or  thermal  radiation  absorbed  and  emitted  by  gases  and  surfaces.  Upward 
longwave  radiative  flux  from  the  ground  is  determined  by  the  surface  emissivity  that  in  turn 
depends  upon  land-use  type,  as  well  as  the  ground  (skin)  temperature.  Shortwave  radiation 
includes  visible  and  surrounding  wavelengths  that  make  up  the  solar  spectrum.  Hence,  the  only 
source  is  the  Sun,  but  processes  include  absorption,  reflection,  and  scattering  in  the  atmosphere 
and  at  surfaces.  For  shortwave  radiation,  the  upward  flux  is  the  reflection  due  to  surface  albedo. 
Within  the  atmosphere  the  radiation  responds  to  model-predicted  cloud  and  water  vapor  dis¬ 
tributions,  as  well  as  specihed  carbon  dioxide,  ozone,  and  (optionally)  trace  gas  concentrations. 
All  the  radiation  schemes  in  WRF  currently  are  column  (one-dimensional)  schemes,  so  each  col¬ 
umn  is  treated  independently,  and  the  fluxes  correspond  to  those  in  inhnite  horizontally  uniform 
planes,  which  is  a  good  approximation  if  the  vertical  thickness  of  the  model  layers  is  much  less 
than  the  horizontal  grid  length.  This  assumption  would  become  less  accurate  at  high  horizontal 
resolution.  Table  8.5  summarizes  the  basic  features  of  the  radiation  schemes  in  the  ARW. 
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8.6.1  Rapid  Radiative  Transfer  Model  (RRTM)  Longwave 

This  RRTM,  which  is  taken  from  MM5,  is  based  on  Mlawer  et  ah  (1997)  and  is  a  spectral-band 
scheme  using  the  correlated-fc  method.  It  uses  pre-set  tables  to  accurately  represent  longwave 
processes  due  to  water  vapor,  ozone,  CO2,  and  trace  gases  (if  present),  as  well  as  accounting  for 
cloud  optical  depth. 


8.6.2  Eta  Geophysical  Fluid  Dynamics  Laboratory  (GFDL)  Long¬ 
wave 

This  longwave  radiation  scheme  is  from  GFDL.  It  follows  the  simplihed  exchange  method  of 
Fels  and  Schwarzkopf  (1975)  and  Schwarzkopf  and  Fels  (1991),  with  calculation  over  spectral 
bands  associated  with  carbon  dioxide,  water  vapor,  and  ozone.  Included  are  Schwarzkopf  and 
Fels  (1985)  transmission  coefficients  for  carbon  dioxide,  a  Roberts  et  ah  (1976)  water  vapor 
continuum,  and  the  effects  of  water  vapor-carbon  dioxide  overlap  and  of  a  Voigt  line-shape 
correction.  The  Rodgers  (1968)  formulation  is  adopted  for  ozone  absorption.  Clouds  are  ran¬ 
domly  overlapped.  This  scheme  is  implemented  to  conduct  comparisons  with  the  operational 
Eta  model. 


8.6.3  Eta  Geophysical  Fluid  Dynamics  Laboratory  (GFDL)  Short¬ 
wave 

This  shortwave  radiation  is  a  GFDL  version  of  the  Lacis  and  Hansen  (1974)  parameterization. 
Effects  of  atmospheric  water  vapor,  ozone  (both  from  Lacis  and  Hansen,  1974),  and  carbon 
dioxide  (Sasamori  et  ah,  1972)  are  employed.  Clouds  are  randomly  overlapped.  Shortwave 
calculations  are  made  using  a  daylight-mean  cosine  solar  zenith  angle  over  the  time  interval 
(given  by  the  radiation  call  frequency). 


8.6.4  MM5  (Dudhia)  Shortwave 

This  scheme  is  base  on  Dudhia  (1989)  and  is  taken  from  MM5.  It  has  a  simple  downward 
integration  of  solar  flux,  accounting  for  clear-air  scattering,  water  vapor  absorption  (Lacis  and 
Hansen,  1974),  and  cloud  albedo  and  absorption.  It  uses  look-up  tables  for  clouds  from  Stephens 
(1978). 

8.6.5  Goddard  Shortwave 

This  scheme  is  based  on  Chou  and  Suarez  (1994).  It  has  a  total  of  11  spectral  bands  and 
considers  diffuse  and  direct  solar  radiation  components  in  a  two-stream  approach  that  accounts 
for  scattered  and  reflected  components.  Ozone  is  considered  with  several  climatological  profiles 
available. 


61 


Table  8.6:  Physics  Interactions.  Columns  correspond  to  model  physical  processes:  radiation 
(Rad),  microphysics  (MP),  cumulus  parameterization  (CP),  planetary  boundary  layer/vertical 
diffusion  (PEL),  and  surface  physics  (Sfc).  Rows  corresponds  to  model  variables  where  i  and  o 
indicate  whether  a  variable  is  input  or  output  (updated)  by  a  physical  process. 


Rad 

MP 

CP 

PEL 

Sfc 

Atmospheric 

Momentum 

i 

io 

State  or 

Pot.  Temp. 

io 

io 

io 

io 

Tendencies 

Water  Vapor 

i 

io 

io 

io 

Cloud 

i 

io 

o 

io 

Precip 

i 

io 

o 

Surface 

Longwave  Up 

i 

o 

Fluxes 

Longwave  Down 

o 

i 

Shortwave  Up 

i 

o 

Shortwave  Down 

o 

i 

Sfc  Convective  Rain 

o 

i 

Sfc  Resolved  Rain 

o 

i 

Heat  Flux 

i 

o 

Moisture  Flux 

i 

o 

Surface  Stress 

i 

o 
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8.7  Physics  Interactions 

While  the  model  physics  parameterizations  are  categorized  in  a  modular  way,  it  should  be 
noted  that  there  are  many  interactions  between  them  via  the  model  state  variables  (potential 
temperature,  moisture,  wind,  etc.)  and  their  tendencies,  and  via  the  surface  fluxes.  Table  8.6 
summarizes  how  the  various  physics  processes  interact  in  the  model.  In  the  table,  i  indicates 
that  the  state  variable  or  flux  is  required  input  for  the  physics  scheme,  and  o  indicates  that  the 
tendency  or  flux  is  a  probable  output  of  the  scheme.  It  can  be  seen  that  all  the  physical  schemes 
interact  in  some  way  with  the  surface  physics  (land-surface  models,  and,  potentially,  coupled 
ocean  models).  The  surface  physics,  while  not  explicitly  producing  tendencies  of  atmospheric 
state  variables,  is  responsible  for  updating  the  land-state  variables. 

Note  also  that,  as  mentioned,  the  microphysics  does  not  output  tendencies,  but  updates  the 
atmospheric  state  at  the  end  of  the  model  time-step.  However,  the  rest  of  the  o’s  in  the  upper 
half  of  the  table  are  representative  of  the  physical  tendencies  of  these  variables  in  the  model. 

The  radiation,  cumulus  parameterization,  and  boundary-layer  schemes  all  output  tendencies, 
but  the  tendencies  are  not  added  utill  later  in  the  solver,  so  from  this  perspective  the  order  of  call 
is  not  important.  Moreover,  these  physics  schemes  do  not  have  to  be  called  at  the  same  frequency 
as  each  other  or  the  model  time  step.  When  lower  frequencies  are  used,  their  tendencies  are  kept 
constant  between  calls.  This  is  typically  done  for  the  radiation  schemes,  which  are  too  expensive 
to  call  every  time,  and  for  the  cumulus  schemes,  for  which  it  is  also  not  necessary.  However, 
the  surface/boundary-layer  schemes  are  normally  called  every  step  in  the  ARW  because  this  is 
likely  to  give  the  best  results. 

The  radiation  is  called  hrst  because  of  the  required  radiative  fluxes  that  are  input  to  the 
land-surface  scheme.  The  land-surface  also  requires  rainfall  from  the  microphysics  and  cumulus 
schemes,  but  that  is  from  the  previous  time-step.  The  boundary-layer  scheme  is  necessarily  after 
the  land-surface  scheme  because  it  requires  the  heat  and  moisture  fluxes. 
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Chapter  9 

Variational  Data  Assimilation 


An  introduction  to  the  basic  ideas  of  variational  data  assimilation  and  the  WRF-Var  system  is 
given  in  this  chapter,  followed  by  a  brief  overview  of  recent  major  improvements  to  WRF-Var. 
This  overview  supplements  the  original  description  of  the  three-dimensional  variational  (3D-Var) 
algorithm  found  in  Barker  et  ah  (2004).  One  of  the  most  important  additions  to  WRF-Var  is  a 
new  utility  gen.be,  used  to  calculate  background  error  covariances  for  a  user’s  own  application; 
it  is  discussed  in  the  latter  half  of  this  chapter.  The  WRF-Var  system  is  evolving  rapidly,  and  a 
future  technical  note  will  accompany  the  general  release  of  the  4D-Var  component  of  WRF-Var. 
That  technical  note  will  include  an  extensive  description  of  the  entire  WRF-Var  system. 

9.1  Introduction 

The  basic  goal  of  any  variational  data  assimilation  system  is  to  produce  an  optimal  estimate 
of  the  true  atmospheric  state  at  analysis  time  through  iterative  solution  of  a  prescribed  cost- 
function  (He  et  ah,  1997): 


J(x)  =  J;,(x)  +  Jo(x)  =  ^(x  -  x^)^B  ^(x  -  x'’)  +  ^(y  -  y°)^(E  -|-  F)  ^(y  -  y°).  (9.1) 

The  variational  problem  can  be  summarized  as  the  iterative  minimization  of  (9.1)  to  hnd  the 
analysis  state  x  that  minimizes  J(x).  This  solution  represents  the  a  posteriori  maximum  like¬ 
lihood  (minimum  variance)  estimate  of  the  true  state  of  the  atmosphere  given  the  two  sources 
of  a  priori  data:  the  hrst  guess  (or  background)  x*^  and  observations  y”  (Lorenc,  1986).  The 
£t  to  individual  data  points  is  weighted  by  estimates  of  their  errors:  B,  E,  and  F  are  the  back¬ 
ground,  observation  (instrumental),  and  representivity  error  covariance  matrices,  respectively. 
Represent ivity  error  is  an  estimate  of  inaccuracies  introduced  in  the  observation  operator  H 
used  to  transform  the  gridded  analysis  x  to  observation  space  y  =  H (x)  for  comparison  against 
observations.  This  error  will  be  resolution  dependent  and  may  also  include  a  contribution  from 
approximations  (e.g.,  linearizations)  in  H. 

As  described  in  Barker  et  ah  (2004),  the  particular  variational  data  assimilation  algorithm 
adopted  in  WRF-Var  is  a  model-space,  incremental  formulation  of  the  variational  problem.  In 
this  approach,  observations,  previous  forecasts,  their  errors,  and  physical  laws  are  combined  to 
produce  analysis  increments  x^  ,  which  are  added  to  the  hrst  guess  x'’  to  provide  an  updated 
analysis. 
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Figure  9.1:  Sketch  showing  the  relationship  between  datasets  (circles),  and  algorithms  (rectan¬ 
gles)  of  the  ARW  system. 


Figure  9.1  illustrates  the  relationship  between  WRF-Var,  the  various  datasets,  and  the  other 
components  of  a  typical  NWP  system  (here  ARW).  The  WRF-Var  assimilation  proceeds  as 
described  in  Barker  et  ah  (2004).  A  number  of  recent  upgrades  to  the  WRF-Var  algorithm  will 
be  described  in  Section  9.2. 

The  three  inputs  to  WRF-Var  are: 

a)  First  guess  —  In  cold-start  mode,  this  is  typically  a  forecast /analysis  from  another 
model  interpolated  to  the  ARW  grid  (and  variables)  via  the  WRF  SI  and  real  programs.  In 
cycling  mode,  the  hrst  guess  is  a  short-range  (typically  1-6  hour)  ARW  forecast. 

b)  Observations  y° —  In  the  current  version  of  WRF-Var,  observations  may  be  supplied 
either  in  a  text  (MM5  3D-Var)  format  or  BUFR  format  (but  not  a  combination  of  the  two). 
An  observation  preprocessor  (3DVAR_OBSPROC)  is  supplied  with  the  code  release  to  perform 
basic  quality  control,  assign  observation  errors,  and  reformat  observations  from  the  MM5  little.r 
text  format  into  3D-Var’s  own  text  format.  Details  can  be  found  in  Barker  et  ah  (2003,  2004). 

c)  Background  error  covariances  B —  used  to  dehne  the  spatial  and  multivariate  response  of 
the  analysis  to  an  observation.  In  variational  systems,  these  covariances  are  typically  calculated 
off-line,  and  signihcant  tuning  is  required  to  optimize  performance  for  a  particular  application 
(e.g.,  Ingleby  (2001);  Wu  et  ah  (2002)).  The  amount  of  work  required  to  do  this  satisfactorily 
is  signihcant,  and  should  not  be  underestimated.  In  order  to  assist  the  user,  WRF  developers 


66 


supply  the  following:  i)  a  default  set  of  statistics  used  for  the  initial  set  up  of  a  domain;  ii)  a 
utility  gen_be  (described  in  Section  9.3)  to  process  ensembles  of  forecasts  into  the  appropriate 
control  variable  space;  and  hi)  diagnostic  routines  to  assess  the  accuracy  of  observation  and 
background  error  statistics.  These  routines  include  both  innovation  vector-based  approaches 
(Hollingsworth  and  Lonnberg,  1986)  and  variational  tuning  approaches  (Desroziers  and  Ivanov, 
2001). 

Following  assimilation  of  all  data,  an  analysis  is  produced  that  must  be  merged  with  the 
existing  lateral  boundary  conditions  x'^'^  (described  in  Barker  et  ah  (2003)).  Note:  In  cycling 
mode,  only  the  wrfbdy  lateral  boundary  condition  hies  (x*’^'^)  output  of  Sl/real  are  used,  and 
not  the  wrfinput  initial  condition  hies  (x^).  In  cold-start  mode,  both  are  required. 


9.2  Improvements  to  the  WRF-Var  Algorithm 

9.2.1  Improved  vertical  interpolation 

The  original  WRF  3D-Var  system  described  in  Barker  et  ah  (2004)  used  height  interpolation 
for  all  observation  operators.  If  an  observation  is  reported  as  a  function  of  pressure,  then  height 
is  approximated  using  the  hydrostatic  relation.  This  step  introduces  an  unnecessary  source 
of  error.  The  new  WRF-Var  system  performs  vertical  interpolation  in  terms  of  the  original 
observed  coordinate,  i.e.,  pressure  or  height. 

9.2.2  Improved  minimization  and  “outer  loop” 

The  default  WRF-Var  cost  function  minimization  uses  a  modihed  version  of  the  limited  memory 
Quasi-Newton  Method  (QNM).  Recently,  an  alternative  Conjugate  Gradient  Method  (CGM)  has 
been  implemented.  Unlike  the  QNM  technique,  the  CGM  method  restricts  3D-Var’s  inner  loop 
to  be  completely  linear.  This  limitation  is  dealt  with  through  the  inclusion  of  an  outer  loop 
in  WRF-Var,  the  purpose  of  which  is  to  iterate  towards  nonlinear  solutions  (e.g.,  observation 
operators,  balance  constraints,  and  the  forecast  itself  in  4D-Var)  using  the  WRF-Var  analysis 
from  the  previous  iteration  as  new  background.  The  outer  loop  is  also  used  as  a  form  of 
variational  quality  control  as  follows:  observations  are  rejected  if  their  0-B  values  are  outside 
a  prescribed  range  (typically  several  times  the  observation  error  standard  deviation).  This 
errormax  test  implicitly  assumes  the  rejected  large  0-B  values  are  due  to  a  bad  observation 
(O)  rather  than  poor  background  (B).  However,  if  it  is  the  background  B  that  is  incorrect  then 
the  system  will  reject  the  most  useful  observations  available  to  the  assimilation  system,  i.e., 
those  in  areas  where  the  hrst-guess  is  poor.  The  outer  loop  alleviates  this  effect  by  allowing 
observations  rejected  in  previous  iterations  to  be  accepted  if  their  new  0-B  falls  within  the 
required  range  in  subsequent  outer  loops.  The  assimilation  of  nearby  observations  in  previous 
iterations  essentially  provides  a  “buddy  check”  to  the  observation  in  question. 

9.2.3  Flexible  choice  of  control  variables 

In  practical  variational  data  assimilation  schemes,  the  background  error  covariance  matrix  B  is 
computed  not  in  model  space  x'  :  u,  u,  T,  g,Ps,  but  in  a  control  variable  space  v  related  to  model 
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space  via  the  control  variable  transform  U,  i.e., 


x'  =  Uv  =  UpU^U,,v.  (9.2) 

The  expansion  U  =  UpU.yU/i  represents  the  various  stages  of  covariance  modeling:  horizontal 
correlations  U/i,  vertical  covariances  U^,,  and  multivariate  covariances  Up. 

The  components  of  v  are  chosen  so  that  their  error  cross-correlations  are  negligible,  thus 
permitting  the  matrix  B  to  be  block-diagonalized.  The  many  varying  applications  (high/low 
resolution,  polar /tropical,  etc.)  of  WRF-Var  require  a  flexible  choice  of  background  error  model. 
This  is  achieved  via  a  namelist  option  “cv.options”  as  dehned  in  Table  9.1. 


cv_options 

2 

(original  MM5) 

3 

(NCEP) 

4 

(Global) 

5 

(Regional) 

Analysis 

Increment 

x' 

Change  of 
Variable 

Up 

i’’ ,x'uKy  ,P'su 

Vertical 

Covariances 

u. 

B  =  EAE^ 

RF 

B  =  EAE^ 

Horizontal 

Correlations 

u,. 

RF 

Spectral 

RF 

Control 

Variables 

V 

v(q  j,m) 

v(h  j,fc) 

v(Z,  n,  m) 

v(z,  j,m) 

Table  9.1:  The  dehnitions  of  the  various  stages  of  the  control  variable  transform  given  by  (9.2)  for 
the  unihed  global/regional  WRF-Var  system.  Indices  {i,j,k)  refer  to  grid-point  space,  index  m 
to  vertical  mode,  and  I,  n  to  global  spectral  mode.  The  variables  are:  u,  v:  velocity  components; 
T:  temperature;  q:  mixing  ratio;  Ps-  surface  pressure;  -0:  streamfunction;  y:  velocity  potential; 
r:  relative  humidity.  The  subscript  u  indicates  an  unbalanced  held.  The  acronym  RF  stands 
for  recursive  hlter. 

Table  9.1  indicates  that  the  only  difference  between  global  (cv_options=4)  and  WRF  regional 
(cv_options=5)  versions  of  the  WRF-Var  control  variable  transform  is  in  the  horizontal  error 
correlations  U/^.  Note  also,  the  only  difference  between  the  old  MM5  background  error  model 
(cv_options=2)  and  WRF  regional  (cv_options=5)  is  in  the  Up  transform.  The  former  imposes  a 
dynamical  balance  constraint  via  an  unbalanced  pressure  control  variable  (Barker  et  ah,  2004), 
whereas  in  the  new  regional  covariance  model,  balance  is  imposed  via  statistical  regression  (see 
Section  9.3  for  details).  This  choice  of  control  variables  is  considered  more  appropriate  for  the 
mass-based  ARW  solver. 
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9.2.4  First  Guess  at  Appropriate  Time  (FGAT) 

A  First  Guess  at  Appropriate  Time  (FGAT)  procedure  has  been  implemented  in  WRF-Var  (Lee 
et  ah,  2004).  The  FGAT  capability  results  in  a  more  accurate  calculation  of  innovation  vectors 
(observation  minus  hrst  guess  difference),  and  hence  a  better  use  of  observations  when  their  valid 
time  differs  from  that  of  the  analysis.  FGAT  is  most  effective  for  the  analysis  of  observations 
from  asynoptic,  moving  platforms  (e.g.,  aircraft  and  satellite  data).  Surface  observations  with 
high  temporal  resolution  also  beneht  from  the  use  of  FGAT. 

9.2.5  Radar  Data  Assimilation 

Numerous  modihcations  have  been  made  in  order  to  assimilate  Doppler  radar  radial  velocity 
and  reflectivity  observations.  Firstly,  vertical  velocity  increments  are  included  in  WRF-Var  via 
the  “Richarson  balance  equation”  that  combines  the  continuity  equation,  adiabatic  thermody¬ 
namic  equation,  and  hydrostatic  relation.  Linear  and  adjoint  codes  of  Richardson’s  equation 
have  been  incorporated  into  WRF-Var.  In  order  to  develop  a  capability  for  Doppler  reflectivity 
assimilation,  we  use  the  total  water  as  a  control  variable,  requiring  a  partitioning  of  the  moisture 
and  water  hydrometeor  increments.  A  warm-rain  parameterization  is  also  included,  which  in¬ 
cludes  condensation  of  water  vapor  into  cloud,  accretion  of  cloud  by  rain,  automatic  conversion 
of  cloud  to  rain,  and  evaporation  of  rain  to  water  vapor.  Finally,  the  observation  operators  for 
Doppler  radial  velocity  and  reflectivity  are  included  in  WRF-Var.  Further  details  and  results  of 
the  radial  velocity  work  can  be  found  in  Xiao  et  ah  (2005).  The  radar  reflectivity  approach  will 
be  described  in  a  future  paper. 

9.2.6  Unified  Regional/Global  3D-Var  Assimilation 

There  are  many  benehts  to  having  a  single  data  assimilation  system  for  regional  and  global 
applications.  The  majority  of  the  code  is  common  to  both  (observation  operators,  minimization, 
background  error  preconditioning,  interpolation,  etc.).  Technically,  the  main  changes  required 
to  extend  the  regional  application  to  global  are  related  to  the  presence  of  a)  the  polar  singularity, 
and  b)  periodic  East-West  boundary  conditions.  Of  course,  there  are  also  scientific  questions 
concerning  the  optimal  mix  of  observations  required  for  global/regional  models,  and  the  choice 
of  control  variables  and  balance  constraints.  A  unihed  global/regional  3D-Var  system  should 
therefore  be  flexible  to  a  variety  of  thinning/quality-control  algorithms  and  also  to  alternative 
formulations  of  the  background  error  covariance  matrix.  This  flexibility  has  been  a  key  design 
requirement  throughout  the  WRF-Var  project. 

The  major  difference  between  regional  and  global  options  in  WRF-Var  is  in  the  dehnition 
of  horizontal  background  error  covariances.  In  regional  mode,  recursive  hlters  (Purser  et  ah, 
2003)  are  used  to  project  observation  information  away  from  the  observation  location.  In  global 
mode,  a  spectral  decomposition  is  applied.  The  main  benehts  of  the  spectral  technique  are  a) 
the  isotropic  and  homogeneous  covariances  that  are  implied  neatly  solve  the  problems  associated 
with  the  pole  (the  pole  is  not  a  special  point  in  spectral  space),  and  b)  horizontal  correlations 
are  dehned  in  terms  of  a  single  function —  the  power  spectrum  (a  function  of  total  wavenumber). 
However,  the  isotropy  of  correlation  dehned  in  spectral  space  is  also  a  weakness —  anisotropies 
need  to  be  dehned  in  an  alternative  manner.  One  solution  to  this  problem  is  to  replace  the 
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spectral  correlations  with  grid-point  correlations  (e.g.,  in  the  Gridpoint  Statistical  Interpolation 
scheme  nnder  development  at  NCEP).  An  alternative  techniqne  is  to  snpplement  the  isotropic 
spectral  correlations  with  an  anisotropic  component  derived  via  grid  transformations,  additional 
control  variables  or  4D-Var.  Research  nsing  the  latter  techniqnes  is  nnderway  using  the  WRF- 
Var  system. 

The  WRF-Var  code  has  been  adapted  to  perform  assimilation  on  a  global,  regular  latitude- 
longitude  grid.  The  major  modifications  are  as  follows. 

a)  Spectral  to  grid-point  transformations  (and  their  adjoints)  have  been  included  in  3D-Var 
to  represent  the  horizontal  component  (U^)  of  the  background  error  covariance  model. 

b)  A  new  global  WRF-Var  registry  was  created  to  accommodate  the  output  of  global  forecast 
models  (currently  only  the  Korean  Meteorological  Administration’s  (KMA)  global  model  has 
been  tested).  The  hnal  analysis  increments  are  written  in  binary  format  and  added  back  to  the 
global  hrst  guess  to  provide  the  global  analysis  in  the  native  model  format. 

c)  Changes  have  been  made  to  allow  for  periodic  boundary  conditions  in  the  Fast-West 
direction. 

d)  A  number  of  minor  changes  have  been  made  to  treat  the  polar  rows  as  special  points  (e.g., 
in  the  grid-point  io  u,v  wind  conversion  in  the  Up  transform  and  the  observation  operators 
for  polar  winds). 


9.3  Background  Error  Covariances 

Forecast  ( “hrst  guess”  or  “background” )  error  covariances  are  a  vital  input  to  variational  data 
assimilation  systems.  They  inhuence  the  analysis  ht  to  observations  and  also  completely  dehne 
the  analysis  response  away  from  observations.  The  latter  impact  is  particularly  important 
in  data-sparse  areas  of  the  globe.  Unlike  ensemble  hlter  data  assimilation  techniques  (e.g., 
the  Ensemble  Adjustment  Kalman  Filter,  the  Ensemble  Transform  Kalman  Filter),  3/4D-Var 
systems  do  not  implicitly  evolve  forecast  error  covariances  in  real-time.  Instead,  climatologic 
statistics  are  usually  estimated  offline.  The  “NMC-method” ,  in  which  forecast  error  covariances 
are  approximated  using  forecast  difference  (e.g.,  T-l-48  minus  T-l-24)  statistics,  is  a  commonly 
used  approach  (Parrish  and  Berber,  1992).  Experiments  at  ECMWF  (Fisher,  2003)  indicate 
superior  statistics  may  be  obtained  using  a  cycling  analysis/forecast  ensemble  prediction  system 
based  on  perturbed  observations/physics. 

Recent  advances  permit  the  use  of  flow-dependent  forecast  error  covariances  in  3D/4D-Var 
through,  for  example,  grid  transformations  (Desroziers,  1997),  anisotropic  recursive  filters  (Wu 
et  ah,  2002;  Purser  et  ah,  2003),  or  observation-space  formulations  of  the  variational  problem 
(Daley  and  Barker,  2001).  Flow-dependence  may  be  enhanced  in  4D-Var  through  the  use  of  the 
forecast  model  to  provide  dynamical  consistency  to  the  evolving  forecast  state  during  4D-Var’s 
time  window  (Rabier  et  ah,  1998).  Still,  the  practical  effort  required  to  specify  and  implement 
flow-dependent  error  covariances  in  3/4D-Var  is  significant. 

The  NMC-method  code  developed  for  MM5  3D-Var  (Barker  et  ah,  2004)  is  nearing  the  end 
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of  its  useful  life.  The  development  of  a  unified  global/regional  WRF-Var  system,  and  its  appli¬ 
cation  to  a  variety  of  models  (e.g.,  ARW,  MM5,  KMA  global  model,  Taiwan’s  Nonhydrostatic 
Forecast  System  [NFS])  has  required  a  new,  efficient,  portable  forecast  background  error  covari¬ 
ance  calculation  code  to  be  written.  There  is  also  a  demand  for  such  a  capability  to  be  available 
and  supported  for  the  wider  3/4D-Var  research  community  for  application  to  their  own  geo¬ 
graphic  areas  of  interest  (the  default  statistics  supplied  with  the  WRF-Var  release  are  designed 
only  as  a  starting  point).  In  this  section,  the  new  gen_be  code  developed  by  NCAR/MMM  to 
generate  forecast  error  statistics  for  use  with  the  WRF-Var  system  is  described. 

The  background  error  covariance  matrix  is  defined  as 

B  =  ^  ~  x'x'^,  (9.3) 

where  the  overbar  denotes  an  average  over  time  and/or  geographical  area.  The  true  background 
error  e  is  not  known  in  reality,  but  is  assumed  to  be  statistically  well-represented  by  a  model  state 
perturbation  x'.  In  the  standard  NMC-method  (Parrish  and  Berber,  1992),  the  perturbation 
x'  is  given  by  the  difference  between  two  forecasts  (e.g.,  24  hour  minus  12  hour)  verifying 
at  the  same  time.  Climatological  estimates  of  background  error  may  then  be  obtained  by 
averaging  such  forecast  differences  over  a  period  of  time  (e.g.,  one  month).  An  alternative 
strategy  proposed  by  (Fisher,  2003)  makes  use  of  ensemble  forecast  output,  defining  the  x' 
vectors  as  ensemble  perturbations  (ensemble  minus  ensemble  mean).  In  either  approach,  the 
end  result  is  an  ensemble  of  model  perturbation  vectors  from  which  estimates  of  background 
error  may  be  derived.  The  new  genJ)e  utility  has  been  designed  to  work  with  either  forecast 
difference,  or  ensemble-based,  perturbations. 

As  described  above,  the  WRF-Var  background  error  covariances  are  specified  not  in  model 
space  x',  but  in  a  control  variable  space  v,  which  is  related  to  the  model  variables  (e.g.,  wind 
components,  temperature,  humidity,  and  surface  pressure)  via  the  control  variable  transform 
dehned  in  (9.2).  Both  (9.2)  and  its  adjoint  are  required  in  WRF-Var.  In  contrast,  the  back¬ 
ground  error  code  performs  the  inverse  control  variable  transform  v  =  in  order 

to  accumulate  statistics  for  each  component  of  the  control  vector  v. 

Using  the  NMC-method,  x'  =  xt2  —  xti  where  T2  and  T1  are  the  forecast  difference  times 
(e.g.,  48h  minus  24h  for  global,  24h  minus  12h  for  regional).  Alternatively,  for  an  ensemble-based 
approach,  x^'  =  x^  —  x,  where  the  overbar  is  an  average  over  ensemble  members  k  =  1,  Ug.  The 
total  number  of  binary  files  produced  by  this  stage  is  uj  x  Ug  where  Uf  is  the  number  of  forecast 
times  used  (e.g.,  for  30  days  with  forecast  every  12  hours,  Uf  =  60).  Using  the  NMC-method, 
Ug  =  1  (1  forecast  difference  per  time).  For  ensemble-based  statistics,  Ug  is  the  number  of 
ensemble  members. 

The  background  error  covariance  generation  code  gen^be  is  designed  to  process  data  from  a 
variety  of  regional/global  models  (e.g.,  ARW,  MM5,  KMA  global  model,  NFS,  etc.),  and  process 
it  in  order  to  provide  error  covariance  statistics  for  use  in  variational  data  assimilation  systems. 
The  initial,  model-dependent  “stage  0”  is  illustrated  in  Fig.  9.2. 

Alternative  models  use  different  grids,  variables,  data  formats,  etc.,  and  so  initial  converters 
are  required  to  transform  model  output  into  a  set  of  standard  perturbation  helds  (and  metadata), 
and  to  output  them  in  a  standard  binary  format  for  further,  model-independent  processing.  The 
standard  grid-point  fields  are  as  follows. 

•  Perturbations:  Streamfunction  -(/'(i,  j,  k),  velocity  potential  x'ih  j)  ^))  temperature  T'{i,  j,  k), 
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BE  Generation:  Stage  0 


Figure  9.2:  Sketch  of  the  role  of  Stage  0  converters  in  transforming  model-specihc  data  (e.g., 
ARW,  KMA  global  model,  etc.)  to  standard  perturbation  helds  and  relevant  metadata  (e.g., 
latitude,  height,  land/sea,  etc.). 

relative  humidity  k),  surface  pressure 

•  Full-helds:  height  z{i,j,k),  latitude  (These  are  required  for  the  production  of 

background  error  statistics  stored  in  terms  of  physics  variables,  rather  than  tied  to  a 
specihc  grid.  This  flexibility  is  included  in  gen_be  through  a  namelist  option  to  dehne  the 
bins  over  which  data  is  averaged  in  a  variety  of  ways  (e.g.,  latitude  height,  grid  points). 
Land-sea  and  orographic  effects  may  also  be  represented  in  this  way. 

In  general,  the  stage.O  converters  are  developed  and  maintained  by  those  supporting  indi¬ 
vidual  models.  Only  the  WRF-to-standard-£elds  converter  gen.be.stageO.wrf  is  maintained  and 
supported  by  the  ARW  effort. 

9.3.1  Removal  of  time- mean 

In  order  to  calculate  covariances  between  helds,  the  average  value  must  hrst  be  removed.  This 
is  performed  in  the  hrst  stage  utility  geri-bestagel. 

9.3.2  Multivariate  Covariances:  Regression  coefficients  and  unbal¬ 
anced  variables 

The  WRF-Var  system  permits  a  variety  of  background  error  covariance  models  to  be  employed, 
as  described  in  Section  9.2.3  above.  The  utility  gen_be  is  used  to  provide  background  error 
statistics  only  for  cv.options  4  and  5. 
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The  second  stage  of  gen^be  (gen_be^stage2)  provides  statistics  for  the  unbalanced  helds  Xu, 
T„,  and  Pgu  used  as  control  variables  in  WRF-Var.  The  unbalanced  control  variables  are  dehned 
as  the  difference  between  full  and  balanced  (or  correlated)  components  of  the  held.  In  this  stage 
of  the  calculation  of  background  errors,  the  balanced  component  of  particular  helds  is  modeled 
via  a  regression  analysis  of  the  held  using  specihed  predictor  helds  (e.g.,  streamfunction;  see 
Wu  et  ah  (2002)  for  further  details).  The  resulting  regression  coefficients  are  output  for  use 
in  WRF-Var’s  Up  transform.  Currently,  three  regression  analyses  are  performed  resulting  in 
three  sets  of  regression  coefficients  (Note:  The  perturbation  notation  has  been  dropped  for  the 
remainder  of  this  chapter  for  clarity.): 

•  Velocity  potential/streamfunction  regression:  Xb  =  ci’] 

•  Temperature/streamfunction  regression:  Th^ki  =  J2k2^ki,k2'4’k2]  and 

•  Surface  pressure/streamfunction  regression:  psb  =  '^kWk'ilJk- 

Data  is  read  from  all  nj  x  rie  hies  and  sorted  into  bins  dehned  via  the  namelist  option 
binjtype.  Regression  coefficients  G{kl,k2)  and  W{k)  are  computed  individually  for  each  bin 
(bin_type=l  is  used  here,  representing  latitudinal  dependence)  in  order  to  allow  representation 
of  diherences  between,  for  example,  polar,  mid-latitude,  and  tropical  dynamical  and  physical 
processes.  In  addition,  the  scalar  coefficient  c  used  to  estimate  velocity  potential  errors  from 
those  of  streamfunction  is  calculated  as  a  function  of  height  to  represent,  for  example,  the 
impact  of  boundary-layer  physics.  Latitudinal/height  smoothing  of  the  resulting  coefficients 
may  be  optionally  performed  to  avoid  artihcial  discontinuities  at  the  edges  of  latitude/height 
boxes. 

Having  computed  regression  coefficients,  the  unbalanced  components  of  the  helds  are  calcu¬ 
lated  asxu  =  X-  ciJ,  =  Tfei  -  Yjk2  Gki,k2i^k2,  and  Psu  =  Ps-Y.k  ^k'tpk-  These  helds  are 
output  for  the  subsequent  calculation  of  the  spatial  covariances  as  described  below. 

9.3.3  Vertical  Covariances:  Eigenvectors/eigenvalues  and  control  vari¬ 
able  projections 

The  third  stage  [genJbestageS)  of  genJbe  calculates  the  statistics  required  for  the  vertical  com¬ 
ponent  of  the  control  variable  transform.  This  calculation  involves  the  projection  of  3D  helds 
on  model-levels  onto  empirical  orthogonal  functions  (EOFs)  of  the  vertical  component  of  back¬ 
ground  error  covariances  (Barker  et  ah,  2004).  For  each  3D  control  variable  {'ip,  Xu,  Tu,  and  r), 
the  vertical  component  of  B,  is  calculated  and  an  eigenvector  decomposition  performed.  The 
resulting  eigenvectors  E  and  eigenvalues  A  are  saved  for  use  in  WRF-Var. 

The  genJ)e  code  calculates  both  domain-averaged  and  local  values  of  the  vertical  component 
of  the  background  error  covariance  matrix.  The  dehnition  of  local  again  depends  on  the  value 
of  the  namelist  variable  bin.type  chosen.  For  example,  for  bin_type=l,  &  kz  x  kz  (where  kz 
is  the  number  of  vertical  levels)  vertical  component  of  B  is  produced  at  every  latitude  (data 
is  averaged  over  time  and  longitude)  for  each  control  variable.  Eigendecomposition  of  the  re¬ 
sulting  climatological  vertical  error  covariances  B  =  EAE^  results  in  both  domain-averaged 
and  local  eigenvectors  E  and  eigenvalues  A.  Both  sets  of  statistics  are  included  in  the  dataset 
supplied  to  WRF-Var,  allowing  the  choice  between  homogeneous  (domain-averaged)  or  local 
(inhomogeneous)  background  error  variances  and  vertical  correlations  to  be  chosen  at  run  time 
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(Barker  et  al.,  2004).  Having  calculated  and  stored  eigenvectors  and  eigenvalues,  the  final  part 
of  gen_be.stage3  is  to  project  the  entire  sequence  of  3D  control  variable  fields  into  EOF  space 
Vv  =  t/-'vp  =  A-V2E^Vp. 

9.3.4  Horizontal  Covariances:  Recursive  filter  lengthscale  (regional), 
or  power  spectra  (global) 

The  last  aspect  of  the  climatological  component  of  background  error  covariance  data  required 
for  WRF-Var  is  the  horizontal  error  correlations,  the  representation  of  which  forms  the  largest 
difference  between  running  WRF-Var  in  regional  and  global  mode.  (It  is  however,  still  a  fairly 
local  change.) 

In  a  global  application  {gen.hestage^-global) ,  power  spectra  are  computed  for  each  of  the  kz 
vertical  modes  of  the  3D  control  variables  if),  Xu,  Tu,  and  r,  and  for  the  2D  control  variable  Psu 
data.  In  contrast,  in  regional  mode,  horizontal  correlations  are  computed  between  grid-points 
of  each  2D  field,  binned  as  a  function  of  distance.  A  Gaussian  curve  is  then  fitted  to  the  data 
as  described  in  Barker  et  al.  (2004)  to  provide  correlation  lengthscales  for  use  in  the  recursive 
filter  algorithm. 
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Appendix  A 

Physical  Constants 


The  following  is  a  list  of  physical  constants  used  in  the  model. 


TT 

= 

3.1415926 

Pi 

k 

= 

0.4 

Von  Karman  constant 

Te 

= 

6.370  X  10®  m 

Radius  of  earth 

9 

= 

9.81  ms“^ 

Acceleration  due  to  gravity 

kle 

= 

7.2921  X  10“®  s-^ 

Angular  rotation  rate  of  the  earth 

ctb 

= 

5.67051  x  10“®  Wm-^K-^ 

Stefan  —  Boltzmann  constant 

Rd 

= 

287  Jkg-iK-i 

Gas  constant  for  dry  air 

Rv 

= 

461.6  Jkg-^K-i 

Gas  constant  for  water  vapor 

Cp 

= 

7  X  Rd/2  J  kg-i  K-i 

Specihc  heat  of  dry  air  at  constant  pressure 

Cu 

= 

Cp-Rd  Jkg-iK-i 

Specihc  heat  of  dry  air  at  constant  volume 

^pv 

= 

4  X  J  kg-i  K-i 

Specihc  heat  of  water  vapor  at  constant  pressure 

^vv 

= 

Cpy  -  Ry  J  kg“^  K“^ 

Specihc  heat  of  water  vapor  at  constant  volume 

^liq 

= 

4190  Jkg-^K-i 

Specihc  heat  capacity  of  water 

^ice 

= 

2106  Jkg-^K-i 

Specihc  heat  capacity  of  ice 

Ly 

= 

2.5x10®  Jkg-i 

Latent  heat  of  vaporization 

Ls 

= 

2.85x10®  Jkg-i 

Latent  heat  of  sublimation 

Lf 

= 

3.50x10®  Jkg-i 

Latent  heat  of  fusion 

Pw 

= 

1.0  X  10®  kg  m“® 

Density  of  liquid  water 
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Appendix  B 
List  of  Symbols 


Symbols  used  in  this  document  are  listed  in  alphabatical  order  in  this  appendix. 


Symbols 

Definition 

a 

A 

B 

c 

Cs 

Ck 

Cr 

Cr 

^  '  max 

Crp 

Cs 

D 

Dnm 

e 

generic  variable 

coefficient  (Chapter  4),  base-state  lapse  rate  constant  (Chapter  5) 
background  error  covariance  matrix 
scalar  coefficient 
speed  of  sound 

a  constant  used  in  TKE  closure 

Courant  number 
maximum  Courant  number 

Courant  number  from  Table  3.1 

activation  Courant  number  in  vertical  velocity  damping 
a  constant  used  in  eddy  viscosity  calculation 
deformation 

deformation  tensor,  where  n,  m  =  1,  2  and  3 

cosine  component  of  the  Coriolis  term  (Chapters  2,  3);  turbulent  kinetic  energy 
(Chapter  4) 

E 

/ 

F 

F 

Fx 

■^cor 

-^1,2 

9 

Gk 

H 

J 

Kdh,dv 

Kh,v 

observation  error  covariance  matrix 
sine  component  of  the  Coriolis  term 
forcing  terms  for  U,  V,  W,  0  and  Qm 
representivity  error  covariance  matrix 

Coriolis  forcing  terms  for  X  =  U,  V,  and  W 

coefficients  for  weighting  functions  in  specihed  boundary  condition 

acceleration  due  to  gravity 

regression  coefficient 

observation  operator 

cost  function 

horizontal  and  vertical  eddy  viscosity  for  gravity  wave  absorbing  layer 
horizontal  and  vertical  eddy  viscosities 
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Symbols  Definition 

Iq  minimum  length  scale  for  dissipation 

lh,v  horizontal  and  vertical  length  scales  for  turbulence 

Icr  critical  length  scale  for  dissipation 

L  latent  heat  of  condensation 

Lx,y  periodicity  length  in  x  and  y 

m  map  scale  factor 

Hs  ratio  of  the  RK3  time  step  to  the  acoustic  time  step 

N  Brunt- Vaisala  frequency 

p  pressure 

p'  perturbation  pressure 

Po  reference  sea-level  pressure 

Ph  hydrostatic  pressure 

Vht,hs  hydrostatic  pressure  at  the  top  and  surface  of  the  model 

Pdht,dhs  dry  hydrostatic  pressure  at  the  top  and  surface  of  the  model 

Ps  surface  pressure 

Pr  Prandtl  number 

q  generic  scalar 

qc,i,r,s  mixing  ratios  for  cloud  water,  ice,  rain  water  and  snow 

qm  generic  mixing  ratios  for  moisture 

qy  mixing  ratio  for  water  vapor 

qys  saturation  mixing  ratio  for  water  vapor 

Qm  generic  coupled  moisture  variable 

r  relative  humidity 

Te  radius  of  earth 

R  remaining  terms  in  equations 

Rd  gas  constant  for  dry  air 

Ry  gas  constant  for  water  vapor 

t  time 

At  a  full  model  time  step 

T  temperature 

To  reference  sea-level  temperature 

u  horizontal  component  of  velocity  in  x-direction 

U  coupled  horizontal  component  of  velocity  in  x-direction  (Chapters  2,  3,  6,  7) 

control  variable  transform  (Chapter  9) 

Uh  horizontal  correlation 

Up  multivariate  covariance 

Uy  vertical  covariance 

V  horizontal  component  of  velocity  in  ^/-direction 

V  three  dimensional  vector  velocity 

V  coupled  horizontal  component  of  velocity  in  ^/-direction 

V  three  dimensional  coupled  vector  velocity 

w  vertical  component  of  velocity 

W  coupled  vertical  component  of  velocity 
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Symbols  Definition 

Wk  regression  coefficient 

^  height 

Zd  depth  of  damping  layer 

ztop  height  of  model  top 

a  inverse  density  of  air 

a'  perturbation  inverse  density  of  air 

a  inverse  density  of  air  for  the  reference  state 

ad  inverse  density  of  dry  air 

ar  local  rotation  angle  between  j/-axis  and  the  meridian 

off-centering  coefficient  for  semi-implicit  acoustic  step 
7  ratio  of  heat  capacities  for  dry  air  at  constant  pressure  and  volume 

7d  divergence  damping  coefficient 

7e  external  mode  damping  coefficient 

7g  damping  coefficient  for  upper  boundary  gravity  wave  absorbing  layer 

7,.  Rayleigh  damping  coefficient 

e  molecular  weight  of  water  over  the  molecular  weight  of  dry  air  (Chapter  4); 

true  background  error  (Chapter  9) 

Tj  terrain-following  hydrostatic-pressure  vertical  coordinate 

fj  contravariant  ‘vertical’  velocity  or  coordinate  velocity 

6  potential  temperature 

6e  equivalent  potential  temperature 

9m  moist  potential  temperature 

0  coupled  potential  temperature 

fi  hydrostatic  pressure  difference  between  surface  and  top  of  the  model 

Ji  reference  state  hydrostatic  pressure  difference  between  surface  and  top  of  the 

model 

Hd  dry  hydrostatic  pressure  difference  between  surface  and  top  of  the  model 

r  acoustic  time  (Chapter  3),  vertical  structure  function  for  Rayleigh  damping 

(Chapter  4) 

Tnm  stress  tensor  (Chapter  4)  where  n.m  =  1,  2  and  3 
Ar  acoustic  time  step 

(f)  geopotential  (Chapters  2,  3,  5);  latitude  (Chapter  9) 

(f)  geopotential  for  reference  state 

perturbation  geopotential 
$  generic  prognostic  variable  (coupled) 

generic  variable  (Chapter  6) 
streamfunction  increment 
x'  velocity  potential  increment 

to  same  as  rj 

coupled  coordinate  velocity 
ffg  angular  rotation  rate  of  the  earth 
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Subscripts/Superscripts  Definition 


Od 

Oh 

Go 

0 

0' 

()** 

0" 


dry 

hydrostatic 

base  state  sea-level  constant 
reference  state 

perturbation  from  reference  state 
full  value  at  a  Runge-Kutta  step 

perturbation  from  Runge-Kutta  step  value  in  acoustic  steps 
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Appendix  C 
Acronyms 


AFWA 

API 

ARPS 

ARW 

BUFR 

CAPE 

CAPS 

CGM 

COAMPS 

COMET 

DTC 

ECMWF 

EOF 

ESMF 

FAA 

FGAT 

FSL 

GFDL 

GFS 

GRIB 

KMA 

LSM 

MRS 

MM5 

MMM 

MRF 

NAM 

NCAR 

NCEP 

NFS 

NMM 

NOAA 

NRL 


Air  Force  Weather  Agency 
Application  Program  Interface 
Advanced  Regional  Prediction  System 
Advanced  Research  WRF 

Binary  Universal  Form  for  the  Representation  of  Meteorological  Data 
Convectively  Available  Potential  Energy 
Center  for  the  Analysis  and  Prediction  of  Storms 
Conjugate  Gradient  Method 

Coupled  Ocean  /  Atmosphere  Mesoscale  Prediction  System 

Cooperative  Program  for  Operational  Meteorology,  Education,  and  Training 

Developmental  Testbed  Center 

The  European  Centre  for  Medium-Range  Weather  Forecasts 

Empirical  Orthogonal  Function 

Earth  System  Modeling  Framework 

Federal  Aviation  Administration 

First  Guess  at  Appropriate  Time 

Forecast  System  Laboratory 

Geophysical  Fluid  Dynamics  Laboratory 

Global  Forecast  System 

Gridded  Binary 

Korean  Meteorological  Administration 
Land  Surface  Model 
Meter  Kilogram  Second 

Pennsylvania  State  /  NGAR  Mesoscale  Model  Version  5 

Mesoscale  and  Microscale  Meteorology  Division 

Medium  Range  Forecast  Model 

North  American  Mesoscale  Model 

National  Genter  for  Atmospheric  Research 

National  Genters  for  Environmental  Prediction 

Non-hydrostatic  Forecast  System  (Gentral  Weather  Bureau  of  Taiwan) 
Nonhydrostatic  Mesoscale  Model 

National  Oceanographic  and  Atmospheric  Administration 
Navy  Research  Laboratory 
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NWP  Numerical  Weather  Prediction 
OSU  Oregon  State  University 

PEL  Planetary  Boundary  Layer 

PPM  Piecewise  Parabolic  Method 

QNM  Quasi  Newton  Method 

RHS  Right  Hand  Side 

RRTM  Rapid  Radiative  Transfer  Model 
RUG  Rapid  Update  Cycle 

SI  (WRF)  Standard  Initialization 

TKE  Turbulent  Kinetic  Energy 

UCAR  University  Corporation  for  Atmospheric  Research 
YSU  Yonsei  University  (Korea) 

VAR  Variational  Assimilation 

WRF  Weather  Research  and  Forecasting  Model 

WSF  WRF  Software  Framework 
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